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Abstract 

In a weakly excitable medium, characterized by a large threshold stimu- 
lus, the free end of an isolated broken plane wave (wave tip) can either ro- 
tate (steadily or unsteadily) around a large excitable core, thereby produc- 
ing a spiral pattern, or retract causing the wave to vanish at boundaries. 
An asymptotic analysis of spiral motion and retraction is carried out in this 
weakly excitable large core regime starting from the free-boundary limit of 
the reaction-diffusion models, valid when the excited region is delimited by 
a thin interface. The wave description is shown to naturally split between 
the tip region and a far region that are smoothly matched on an intermediate 
scale. This separation allows us to rigorously derive an equation of motion for 
the wave tip, with the large scale motion of the spiral wavefront slaved to the 
tip. This kinematic description provides both a physical picture and exact 
predictions for a wide range of wave behavior, including: (i) steady rotation 
(frequency and core radius) , (ii) exact treatment of the meandering instability 
in the free-boundary limit with the prediction that the frequency of unstable 
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motion is half the primary steady frequency (hi) drift under external actions 
(external held with application to axisymmetric scroll ring motion in three- 
dimensions, and spatial or /and time-dependent variation of excitability), and 
(iv) the dynamics of multi-armed spiral waves with the new prediction that 
steadily rotating waves with two or more arms are linearly unstable. Numer- 
ical simulations of FitzHug-Nagumo kinetics are used to test several aspects 
of our results. In addition, we discuss the semi-quantitative extension of this 
theory to finite cores and pinpoint mathematical subtleties related to the thin 

interface limit of singly diffusive reaction-diffusion models. 
PACS numbers : 47.20 Hw, 82.40.Bj, 87.22. As 
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I. INTRODUCTION 



Spiral waves are characteristic structures of excitable media that have been observed 
in systems as different as catalytic surface oxidation ||, the Belousov-Zhabotinsky chem- 
ical reaction @-[7| aggregating colonies of slime mold ||, and heart tissue where they are 
suspected to play an essential role in cardiac arrhythmia and fibrillation ||. Spiral waves 
are prone to a variety of instabilities, the best studied of which is meander, and they can 
be made to drift and be controlled in diverse ways, for instance by varying the medium 
excitability in space or/and time, or by adding an external field. 

Much of the observed experimental phenomenology has been reproduced by using simpli- 
fied two-variable activator-controller types of description, like the classic FitzHugh-Nagumo 



(FN) model []10| and mild variations of it. Extensive surveys |TT| , p!2| of the possible types 
of wave motion in such models have been performed in a reduced parameter space where 
the only two parameters left to vary are the medium excitability A, defined in section II 
in such a way that the isolated pulse speed is proportional to A for weak excitability, and 
the ratio e between the time scale of the activator and controller kinetics, which controls 
the abruptness of the wave front (i.e. the thickness of the interface delimiting the excited 
region). Different regimes have been identified that are summarized in Fig. [I] for simple 
FN kinetics [|T^]. In the whole region above the propagation boundary {dP), the medium 
excitability is too weak for any plane wave to propagate persistently. In the narrower region 
comprised between dP and the rotor boundary (dR), the medium excitability is sufficient 
for a plane wave to propagate but not for a spiral wave to form. In this region, the end of 
a broken wave-front, referred hereafter as the 'wave tip', simply retracts steadily (Fig. ^|a), 
such that this finger-shaped wave must shrink in length and eventually vanish at boundaries 
in a finite system. In the even narrower region comprised between dR and the meander 
boundary (<9M), the excitability is sufficient for large core spirals to form and the wave tip 
now rotates steadily at a frequency u\ around a circular core of radius i?o- Right on the dR 
boundary, a half plane wave, referred hereafter as the "critical finger", propagates without 
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changing its shape. It can be equivalently interpreted as a retracting finger with vanishing 
retracting velocity or as a spiral wave of infinite core radius. As one keeps increasing the 
excitability, the radius of the spiral core decreases and below the dM boundary the spiral 
tip traces a classic "flower-like" meander pattern (Fig. |2|b). It has been shown that meander 
originates from a supercritical Hopf bifurcation at dM which adds a second frequency uo 2 
to the basic spiral rotation |14]J15l| . The meander patterns exhibit first inward petals as 
uji < uj\. Outward petals appear as u 2 becomes greater than oj\. Further away from <9M, 
the spiral tip motion becomes more complex and possibly chaotic passed the still poorly 
characterized dC boundary ]l2j (not shown in Fig. |l|). Given that a Hopf bifurcation takes 
place on dM, symmetry arguments fix its resonant coupling to the translation modes when 
uj 2 = u>i and thus determine the bifurcation structure of the tip motion near the codimension 



2 point lu 2 = 0J\ on dM [16,17 



In contrast to this rather detailed knowledge, the precise mechanisms that govern spi- 
ral formation and motion remain less well understood from both physical and predictive 
viewpoints. A simple picture remains missing to answer even basic questions, such as why 
does meander occur and why is this instability oscillatory (i.e. a Hopf bifurcation) beyond 
numerical observation ? From a predictive viewpoint, we still lack a quantitative analytical 
understanding of what controls the dM boundary or the frequency ratio 102/ oj% in the pa- 
rameter space of react ion- diffusion models. Similar uncertainties are to a large extent also 
present for other phenomena like spiral drift under external action. 

A kinematical model of spiral dynamics, aimed at the weakly excitable large core limit, 



has been proposed some years ago on a purely phenomenological basis |18| . It has been 
helpful to rationalize experimental facts but it has not been derived from the underlying 
react ion- diffusion equations. Thus it remains limited in its predictions, e.g. it falls short of 
predicting the dM boundary and the ratio u^/^i- Moreover, at a more conceptual level, 
the general validity of the boundary condition assumed for the free end of the wave-front in 
this kinematic theory has remains somewhat unclear. 

The first goal of this paper is to present a rigorous asymptotic derivation of a kinematic 



theory of spiral wave motion in the weakly excitable and free-boundary limit (lower left hand 
corner of Fig. 1) onto which we focus. As we shall see, the structure of this theory differs 
from the one proposed phenomenologically in Ref . . The second goal of this paper is to 
demonstrate, through selected applications of this theory, that it is able to provide a physical 
and quantitative understanding of a wide range of wave phenomena such as, meander, drift 
under various external actions considered in previous studies ||l9|-p5[], and multi-arm spiral 
wave motion. Highlights of our results include an asymptotically exact treatment of the 
meander instability for e«l, which gives the precise location of the DM boundary and 
shows that the instability arises from a supercritical Hopf bifurcation with io-iju)\ = 1/2 in 
this limit, the finding that multi-arm spiral waves with two or more arms are always linearly 
unstable, in contrast to a previous numerical study and predictions of the spiral drift 



speed and drift angle in an external field. These results are generally found to be in good 
quantitative agreement with our simulations of FN kinetics. 

The starting point of our analysis is the standard free-boundary limit of reaction-diffusion 
models [E7[ described in section II, which is valid when the excited region is surrounded by 



a thin interface of width e«l. In this limit, the fast activator variable is eliminated in 
favor of an eikonal equation that gives the normal velocity of this boundary. This velocity 
generally depends on the local radius curvature of this interface, assumed large compared 
to e, as well as the local value of the slow controller variable at the interface. 

This free-boundary problem is non-trivial to solve because it requires to treat both the 
dynamics of the wave- front, which is the part of the boundary where the excited region 
propagates into the recovery region of the medium, and the wave-back where the reverse 
process occurs. Far from the tip, the front and back behave essentially identically, such that a 
'single front' description is rigorously possible. In the tip region, however, the front and back 
must be matched at the tip (i.e. point of zero normal velocity along the boundary), which 
is a difficult task. For this reason a single-front description with a somewhat arbitrary tip 
condition was first used historically to relate the steady rotation frequency and core radius 
of spirals f28| . The kinematic theory of Ref. [|l^] is an attempt to extend this picture to 



an unsteady situation. Subsequent solutions of the complete free-boundary problem, with a 
rigorous matching of front and back that provides a unique and independent determination 
of the spiral frequency and wavelength, focused on two limits. One of these limits (see |29| 
and earlier references therein) is obtained mathematically by assuming e< 1 while keeping 
A fixed of order unity, which corresponds physically to a highly excitable medium. The 
wavelength and frequency obey in this limit certain scaling laws with e first proposed by 



Fife The wave- front and wave-back, however, are matched onto singular core solutions 
(of size e with only activator diffusion) that have later been shown to be generically unstable 
; this result actually seems to concord with the numerical observation of complex meander 



(and thus unstable motion) in this limit | 12| . Thus these solutions do not provide a proper 
starting point for a kinematic theory aimed to describe the onset of meander. A better 
starting point, on which we focus here, is the second limit originating from Ref. [B3| where 



one constructs smooth core solutions to the free-boundary problem. It was shown in |33 
that this can in fact only be consistently done for a weakly excitable medium when the 
radius of curvature of the boundary at the tip remains much larger than the front and back 
interface width, i.e. by assuming simultaneously e <C 1 and A ~ e 1 / 3 <C 1. This allowed 
a rigorous derivation of the line OR in the weak excitability limit [[3^] in good quantitative 
agreement with numerical simulations of FitzHugh-Nagumo kinetics (lower left hand corner 
of Fig. 1), as well as a semi-analytic derivation of the selected core radius/frequency of spiral 



waves and retracting wave speed in the same limit |34| . 

The present kinematic theory is derived by first refining analytically the description of 
steady retraction and rotation in this weakly excitable limit (section IV), and then extending 
it to an unsteady regime (section V) for the non-trivial case of self interacting spirals, i.e. 
where the wave tip motion is influenced by the average controller concentration left by the 
previous passage of the wave- front. This allows us to derive an equation of motion for the 
wave tip that is then used to analyze meander in a linear and nonlinear regime (section VI). 
Results of these two sections have been summarized in a previous short publication . 
Further applications of the kinematic theory are then contained in subsequent sections that 
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include spiral drift under various external actions without self-interaction (section VII) and 
interacting multi-arm spiral waves (section VIII). Finally, corrections to the large core results 
are discussed in Section ITX| and several points are further analyzed in four appendices. In 



addition, for clarity of exposition, we have found best to first give a simple physical picture of 
the kinematic theory and summarize the main results of its application in section III. This 
section is purposely aimed to discuss this theory in terms of experimentally measurable 
quantities as well as to provide a guide to the rest of the paper. 

II. REACTION-DIFFUSION MODEL AND FREE-BOUNDARY LIMIT 

We consider the classic activator (u) controller (v ) two- variable reaction-diffusion model 



of excitable media 10 



d t u = D u V 2 u + f(u,v)/r u (1) 
8 t v = D v V 2 v + g(u,v)/r v (2) 

with a linearly stable rest state (uq,Vq). We focus in this paper on the singly diffusive 
case D v = 0, although we shall also briefly consider the slow controller diffusion limit 



7 = D v /D u <C 1 in section [VI]. The u-nullcline (f(u, v) = 0) is assumed to have the standard 
S-shape in the (u, v) plane. A simple choice of FN kinetics that we use for the numerical 
simulations is f(u, v) = 3u — u 3 — v, g(u, v) — u — S with the rest state u = 6, v = 35 — 5 3 . 
It is convenient to rewrite Eqs. ([l|) and (fj) in a standard dimensionless form by measuring 
time and length in units of t v and (D u t 2 /t u Y^ 2 , respectively, which yields for the singly 
diffusive case 

d t u = eV 2 u + f(u,v)/e (3) 
d t v = g(u,v) (4) 

where e = t u /t v . We study this dimensionless form of the equations in the rest of this paper, 
except in the next section where we summarize the essential ingredients of the kinematic 



theory in dimensional units. For small e, the excited region (u ~ \/3 with the previous 
choice) of a propagating wave is separated by a sharp boundary from the unexcited or 
recovering medium {u — Vo with the previous choice). The wave description can thus be 
reduced to determining the motion of its boundary (i.e. a free boundary problem) P J28| , |32 



(5) 

d t v = g(u ± (v),v) in V ± , (6) 

where c n is the normal velocity of the interface separating the excited and recovery regions 
of the medium denoted by (T> + ) and (T>~~), respectively, k is the local curvature of this 
interface, and M ± (f ) denotes the right-most (+) and left-most (— ) branch of the u-nullcline 
(f(u,v) = 0). The function c(v) is entirely determined by Eq. ([3]) with v fixed. We measure 
the excitability of the medium, i.e. the threshold stimulus necessary to cause a response, by 
the parameter A = v s — Vq, where v s is the stall value of v at which c(v s ) = 0. The isolated 
pulse speed Cq = c(vq) is then a monotonously increasing function of A with cq = aA for 
A -C 1 (q = 1/V2 for our numerical choice). For values of of v near v s , Eq. (^|) can be 
simplified even further to 

d t v = l/r e in P+ (7) 
d t v = -(v -v q )/t r in V~ (8) 

where the activator time scale r e = l/g(u + (v s ), v s ) controls the pulse duration and the 
recovery time 

TR = (d u gdj-d v gd u f) lu=u ' {v ^ (9) 

is the time scale over which the controller variable returns to its rest state after an excitation 
(for our numerical choice, r e = 1/(2^/3), t r = 6). 

III. PHYSICAL PICTURE OF KINEMATIC THEORY AND MAIN RESULTS 
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A. Retraction and rotation 



In a typical chemical or biological excitable media, many parameters (chemical or ionic 
concentrations, temperature, light etc.), control the excitability of the medium. However, 
independently of the complexity of the medium, it is generally possible to construct a single 
dimensionless parameter P3"||£l| 

B = ^ = ^> <"» 

which determines whether the tip rotates, and thus forms a spiral wave, or whether it simply 
retracts. This parameter is expressed as the ratio of two lengthscales that characterize the 
tip region of a broken plane wave that is shown schematically in Fig. ^. The first is the 
radius of curvature R t i P of the wave boundary at the tip. In the limit (e <C 1) where this 
boundary is thin, R tip is obtained by applying the eikonal equation at the tip, which yields 
c n = = c — D u / Rap, and thus R tiv = D u /c where c is the plane wave speed. The 
second lengthscale W is the constant width of the excited region away from the tip. As 
argued in |33| , |34|1 , the wave boundary in the tip region can only be smooth if Ru p ~ W, 
such that the OR boundary must correspond to a fixed value of B = B c of order unity; an 
explicit calculation (Refs. p3| , |34|| and section IV.A here) yields B c = 0.535. For B > B c , the 
excitability of the medium is not sufficient to overcome diffusion in the most highly curved 
part of the tip region, which retracts. In contrast for B < B c , the excitability is sufficient 
to overcome retraction, and the increase of c n away from the tip induces rotation. 

Increasing (decreasing) B corresponds to decreasing (increasing) the excitability of the 
medium while moving normal to the dR boundary in the multi-dimensional parameter space 
that characterizes a given excitable medium. B is therefore the most natural parameter to 
characterize this excitability close to this boundary, and is used throughout this paper. From 
an experimental standpoint, Cq and W are in principle measurable quantities and D u can be 
either measured or estimated, such that one could attempt to quantify B directly. Of course, 
in practice, the definition of W, and thus B, becomes less precise when e is not small and it 
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is simpler to use an experimental control parameter P ex (such as a concentration) that can 
be varied to cross the dR boundary. Therefore we will briefly describe in subsection [111 C 



a more direct way to obtain a quantitative relationship between B and P ex close to the dR 
boundary without the need to use Eq. (JlQ|) . 

B. Rigid wave tip and slaved wave-front 

Close to the dR boundary, we will show that a spiral wave in its tip region, behaves as 
a 'rigid body' whose motion can be characterized by giving two instantaneous quantities: 
the tip tangential velocity Ct(t), and its rotation rate u>i(t) = Ct(t)/Ri(t) where Ri(t) is the 
radius of curvature of the tip trajectory as depicted in Fig. [|. Two key ingredients make 
this kinematic description possible. The first is that, near dR, the wave shape is close to 
a critical finger (i.e. the broken plane wave that simply translates without retraction or 



rotation for B = B c ) on a length scale £ ~ D u /(coyl — c t /c ) (as explained in section IV. C) 
that is large compared to the scale of the tip itself Ru p = D u /cq. The second is that the 
relations that govern q and Ri are established on time scales that are both much shorter 
than the steady rotation period, T = 2tc/uji, as discussed at the start of section V and in 
quantitative details in appendix B. This separation of time scales, which makes our adiabatic 
treatment possible, becomes exact in the large core limit B — > B c , but, importantly, it does 
not depend on e being small. Therefore the present kinematic description should rigorously 
extend beyond the free-boundary limit but we shall assume that e C 1 in this paper to 
compute c t . 



A key difference between this description and the one proposed in [18| should already 
be apparent. Namely here, the dynamics is driven entirely by the rigid tip region (which is 
just a point on the core scale). In contrast, in the kinematic model of the tip motion 
is determined by the wave-front dynamics via a boundary condition imposed at the tip. In 
the present context, the dynamics of the spiral wave-front outside the tip region need not be 
invoked to calculate the tip motion. 
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C. Tangential velocity 



The tip tangential velocity is determined by the controller concentration (equivalently the 
spatial variation of excitability in which the wave-front propagates) in the tip neighborhood. 
To compute this velocity we exploit the fact that, close to the dR boundary, the equations 
of motion can be linearized around the critical finger. This allows us to obtain a solvability 
condition (i.e. a general condition for the existence of a solution to these linearized equations) 
that uniquely relates the tangential velocity c t ({v }) to an arbitrary spatial distribution of 
v; this distribution is only constrained to deviate slightly from the rest state vq in order for 
the wave tip to remain close to the critical finger. This solvability condition is first used in 
section IV.B in the simplest steady-state situation where the tip propagates into a uniform 
controller concentration. The result of interest is 

q/c = 1 + (B - B c )/K (11) 

where K ~ 0.63 is a numerical constant. This result implies that on the weak excitability side 
of the dR boundary (B > B c ), steady retracting waves form with q/co = \Jc^ + Cq/cq > 1, 
where c r is the tip retracting speed. On the other side (B < B c ), spiral forms with q/co < 1 
and a second relation discussed in the next subsection is needed in this case to determine 
the rotation rate. The calculation of the tangential velocity is extended to the case of steady 
self-interacting spirals that propagate in a not fully recovered medium in section IV. D, to 
unsteady self-interacting spirals in section V.B, and to spirals in an external field in section 
VII. B. In the latter case, c t ({v}, E\\, E±) depends both on the controller concentration and 
the components En and E± of the external field respectively parallel and orthogonal to q. 

A relation between B and an arbitrary experimental control parameter P ex can be ob- 
tained close to the dR boundary by simply measuring the slope S of the curve c t /c vs P ex , 
which should be the same on both sides of dR and presumably simpler to obtain on the 
retracting side. It then follows at once from Eq. (|TT|) that 

B-B c = KS(P ex -P ex , c ) (12) 
11 



close to dR, where P ex , c is the value of P ex where dR is crossed. This relationship can be 
used to relate quantitatively the results of the rest of this paper to experiments, keeping in 
mind that these results are only accurate asymptotically close to dR and for small e. 



D. Rotation rate 



The motion of the wave tip region, although rigid, must generally be consistent with 
the motion of the rest of the wave-front away from the tip. On the spiral side of dR, the 
tip region must necessarily rotate to accommodate the fact that its tip end translates at a 
slower speed than the plane waves radiated outward from the core (on the other side of dR, 
Ct > Co simply implies retraction of the tip). In section IV. C, we show that the tip and the 
far regions can be matched on the gently curved intermediate scale i, yielding a rotation 
rate c t /Ri, with 



R, 



Du 
c 



3/2 



(13) 



[1 - q/c )_ 

where b ~ 2.946 is a constant that is obtained by matching the curved tip and far regions. 



It should be noted that this constant differs from the constant b' |36[ obtained by arbitrarily 
imposing a radial departure of the wavefront from the steady-state circular core trajectory 
as in Ref. |[L8|| . The 3/2 exponent, however, is the same both here and in Ref. |l8j since it 
does not depend on details of the matching on the tip scale. 

Eq. ([0^) holds both for steady rotation (in the context of which it is derived in section 
IV. C) and for unsteady rotation owing to the aforementioned adiabatic approximation. For 
steady rotation, the core radius R is simply obtained by substituting the expression for c t 
from Eq. ([□]) into Eq. ([IB]) , which yields, 



R(\ 



Du 
c 



bK 



3/2 



(B c -B) 

The generalization to self- interacting spirals is given in section IV. D. 



(14) 
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E. Parameterization of the wave tip trajectory 



Knowing how to compute Ct and Ri gives in principle a complete kinematic theory of 
the wave tip motion, since this uniquely predicts the Euclidean trajectory of the tip in time. 
However to characterize analytically the tip dynamics in unsteady situations (such as drift, 
meander, etc) it is convenient to measure the instantaneous tip position by the standard 
polar coordinates (r, 0) with respect to a fixed origin at the center of steady rotation (Fig. ^|), 
and to relate the tip motion in these coordinates to q and Ri. This part of our analysis is 
carried out in section V.A and yields a simple forced harmonic oscillator equation 

^ +ul5r = u^SRi (15) 

where 5r(t) = r(t) — Ro is the radial displacement of the wave tip from its radius Ro of 
steady rotation and 5Ri(t) = Ri(t) — Ro. Eq. fllBD is valid for a small radial displacement 
(\Sr\/Ro <C 1 ) and is accompanied by an independent equation for the angular displacement 
ijj(t) = 8(t) — uj\t from steady rotation. For a small radial displacement, however, the two 
equations are not coupled such that Sr can be computed independently. Without forcing, the 
solution of Eq. (|15"D, namely an harmonic motion at frequency u\, is a simple superposition 
of the two translation modes: it gives the tip displacement of a steady spiral which is slightly 
translated with respect to the reference unperturbed spiral. 



F. Main results 

In summary, the application of the present kinematic theory contains three steps: (i) 
using a solvability condition to calculate q in terms of the local controller concentration, 
external field, etc., with a resulting expression that depends on the situation considered, (ii) 
using Eq. (|I3"D to express Ri in terms of c t , and (iii) solving Eq. (|T5D to obtain the radial 
displacement of the tip for a given forcing, which also obviously depends on the situation 
considered. We now summarize the result of this procedure for the selected applications 
examined in this paper (in a different order than in subsequent sections). 
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The simplest example (section VILA) is to compute the the tip motion induced by a small 
periodic spatially uniform variation of excitability B(t) = B + 5Bsm(uJit + <fi). Following 
the above steps, and using the fact that the perturbation is small, we obtain at once 

— y + ulSr = cul(dR /dB)5Bsm(cu 1 t + 4>) (16) 

where the function Ro(B) is defined by Eq. flU]) , and dR /dB is to be evaluated at B = B . 
This resonant forced harmonic oscillator equation has a growing sinusoidal solution with an 
amplitude that increases linearly in time, and which thus corresponds to a spiral drift at a 
speed Cd = Ui(dRo/dB)5B/2, or c</ = iO\ dRo/dP ex 5P ex /2 in an experiment. The action of 
an electric field is considered in section [VII B| and produces a similar type of periodic forcing 



that leads to spiral drift. In agreement with previous studies [ p0| , pi| ,p^-p5| the spiral is 
found to drift at an angle with the external field. This result also determines the curvature- 
induced motion of a scroll filament. There the main prediction is that rings expand in the 
large core limit (i.e. the filament tension is negative) in agreement with previous numerical 
observations in this limit (see [37| and earlier references therein). 



In the case of meander the tip tangential velocity and hence the forcing on the right- 
hand-side of Eq. flTSP depends on the radial displacement 5r(t) — 8r(t — T ) of the tip after 
one complete rotation (Fig. |||) due to the self-interaction of the wave-front with its own 
recovery tail. If this displacement is positive, the average controller concentration will be 
slightly more elevated in the tip region (i.e. the medium will be slightly less excitable in this 
region) than if it is negative, which then affects c t ({v}) and thus Ri and the forcing of the 
tip. This effect leads to a differential equation with delay of the form 

^| +ulq = ulmF{q{t)-q{t-T,)) (17) 

where we have defined the dimensionless radial displacement q(t) = Sr(t)/Ru p = co8r(t) / D u , 
the parameter m = 3B c (dRo/dB)e~ T °^ TR , and F is a tanh-shaped function which we compute 
in section V.B. The saturation of F at large radial displacement is due to the fact that the 
controller concentration only varies appreciably on the scale of Ru p - A linear stability 
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analysis of this equation in section VIA yields that the onset of meander occurs when 
m exceeds a threshold 3/[8F'(0)] that depends in a singular way on the diffusivity ratio 
D v /D u and e. Namely, the function F is non-analytic at in the pure sharp boundary 
limit (^],§D of the singly diffusive model, which sheds some light on difficulties that were 
previously encountered when attempting to perform a linear stability analysis in this limit 
fl38fl . However, for a finite interface width y/D u r u , small compared to the spiral tip radius 
D u /co, the slope at the origin is finite, with F'(0) ~ —\ii(co\Jt u /D u ), such that there is a 
finite meander threshold that will be typically of order unity in experiments or simulations. 
In addition, this analysis predicts that U2/0J1 = 1/2 at onset in the large core limit and a 
simple physical interpretation of both the existence of a threshold and oscillatory motion is 
given at the end of section VI. A. Slightly away from the large core limit, the discussion in 
section IX leads to a modified differential equation with delay that shows that 002/ 'oj\ increases 
above 1/2 as the core radius Ro is decreased, in semi-quantitative agreement with numerical 
simulations. This actually provides a simple picture of the onset of quasiperiodicity (i.e. 
how 002/uJi becomes irrational) as one moves away from the large core limit. 

Finally, for spirals with N arms, we obtain a system of N coupled differential equations 
with delay and with the interaction between the arms controlled by the parameter mjy = 
3B c (dRo/dB)e~ T °^ NrR \ An exact linear stability analysis shows that, unlike for meander, 
there is no finite threshold of instability. Moreover, this instability develops on a time scale 
proportional to l/m\ for N = 2 and l/rajy for N > 2, such that the time necessary to 
observe it grows exponentially as the dR boundary is approached. 

IV. STEADY STATES 

We start by analyzing steady wave patterns of the free boundary problem (|^-§|). As 
shown in ref . |32||^|] , when excitability is decreased, the spiral core radius and spiral period 
diverge on the line dR with A = A c (e) which marks the lower excitability limit of spiral wave 
propagation in the (e, A) plane. As described in section I, on the line dR, spirals degenerate 
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into critical fingers that translate at cq, the plane wave speed. For A < A c (e), the steady 
waves are retracting fingers. Laws for the tip retraction speed and spiral tip divergence were 



obtained in ref. |34]] from numerical computations in the neighborhood of the line A c (e) 
for e <C 1. Here, we begin by recalling the result of |33| about the line of existence of 
critical fingers. We then proceed and study steady patterns in the neighborhood of A c (e) 
by perturbation around the critical fingers. On the retracting wave side, we determine the 
tangential speed of the tip as a function of A — A c (e) by a solvability condition [39|. This 
is the simplest example of the method that we will use in more complicated situations to 
determine the tip tangential speed. For A > A c (e), the critical finger winds up around its 
tip and becomes a steady spiral rotating around a circular core Rq at a constant tangential 
speed q. We first consider the case where Rq is large enough so that the spiral front interface 
can be assumed to propagate in the medium rest state (i.e. the disturbance of the medium 
induced by the tip previous passage can be neglected). We obtain analytically the divergence 
of the spiral radius by adding to the previous determination of c t , an analysis of the Burton- 
Cabrera-Frank (BCF) equation in the large radius limit using matched asymptotics, 



thus confirming the laws obtained in ref. |34j| . Finally, we determine the modification of the 
steady spiral parameters induced by the perturbation of the medium characteristics due to 
previous passages of the spiral. 



A. The line A c (e) of critical fingers 

We now examine the critical fingers that propagate in a shape preserving way at the 
pulse speed c on the boundary A c (e) in the parameter space (A,e). For a small medium 
excitability, the scaling of the line A c (e) is easily determined by comparing two lengthscales 
p3| as reviewed in section III. A. Firstly, the condition that the normal velocity vanishes at 
the wave tip requires that the tip radius of curvature is equal to e/co- This gives the order of 
magnitude of the distance between the wave front and back interfaces. Secondly, the front 
and back interfaces should move at the same velocity. The value of the controller field v 
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should therefore increase from v o = v s — A on the front interface to v s + A on the back 
interface in a time (e/c )/c . This gives r e e/cl ~ A and, remembering that c = «A, the 
scaling A c (e) ~ (e/a 2 r e ) 1/3 . 

A more detailed analysis is required to determine the constant in this asymptotic relation 
and the critical finger shape that we will subsequently need. We follow ref. and search 



for a steady state finger shape translating at cq, the isolated pulse speed. It is convenient 
to work in the frame of the finger with the origin at the finger tip (see Fig.|]) and to use as 
length unit e/ Co, the finger tip radius. On the front interface Yf(x), the value of the controller 
field is equal to the rest state value vq. At point x on the back interface, the controller field 
value has increased to v + e(Yf(x) — Y b {x))/c^r e from Eq. (^). Eqs.(||) therefore become : 
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+ [1-B [Y f (x) - Y b {x)]\ 



1 + 



,dY, 



b^2 



dx 



3/2 



(18) 
(19) 



where Yf(x) denotes the front interface of the finger and Y b (x) its back interface. These equa- 
tions depends on the single parameter B = e/(a 2 r e A 3 ) [[y]]. The searched solutions should 
satisfy at the tip the boundary conditions Yf(0) = Y b (0) = 0, dYf/dx(0) = —dY b /dx(0) = 
+oo, and asymptotically dYf/dx(+oo) = dY b /dx(+oo) = 0. 

The solution of fll8D does not require any supplementary condition. At x — 0, it tends to 
zero as Yf(x) ~ \[2x in agreement with the chosen tip radius of length unity. At x = +oo, 
it diverges logarithmically Yf(x) ~ 21n(x). In fact, Yf(x) can be obtained analytically, 

2 



x = 2 arctan(t> ) + 



v-l 



71 



Y) = ln 



v 2 + l 



(20) 



_(v-l) 2 _ 
with 1 < v < +oo. 

On the contrary, Eq.(|T9|) can be solved with the appropriate boundary conditions only 
for a particular value of the parameter B. The two boundary conditions at x = entirely 



determines the solution of ([19]) once the front interface is determined. The solution Y b {x) 
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should approach [Yf(x) — 2/B c ] as x — > +00 to satisfy the boundary condition at infinity. 
A linearization of (|19| ) around this asymptotic behavior gives a convergent mode and a 
divergent mode growing like ex^{\fB~ c x). So, the solution obeys the right boundary condition 
at x = +00 only for the special values of B which cancel the prefactor of the diverging mode. 
This is numerically found to happen for B c = 0.5353 ■ • • which defines the line of existence 
A c (e) of the critical fingers in the (e, A) plane. In the following, we refer to the solution of 
Eq. (|18|) and ([19|) with B = B c as the " critical finger shape" . It is plotted in Fig .f|. 

Remark. On can note that at the level of Eq. ( |l8l) and (|l~9| ) the interface is continuous as 
well as its first two derivatives. However, its third derivative is discontinuous at the finger tip 
(x — 0,y — 0) since one has Yf(x) = \^2x+x/3 + - ■ ■ while Y^x) = — \^2x+x(l — 2B c )/3 + - ■ ■ 
(and B c 7^ 0). This weak non-analyticity can be cured by introducing a small boundary layer 
near the tip as discussed in Appendix [A|. 



B. Retracting fingers 



We consider a medium characterized by a parameter B = e/(a 2 r e A 3 ) higher than B c , 
that is, not excitable enough to allow for the existence of spirals. We look for steady state 
shapes propagating at c t . We use as before e/co as unit length where Co = a A is the velocity 
of planar front in the considered medium. Eqs.(|5|) determining the front yj and back y b 
interfaces become, 

3/2 
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(21) 
(22) 



The solutions should satisfy at the tip the same boundary conditions as the critical fingers, 
= 2/6(0) = 0, dyf/dx(0) = —dyb/dx(0) = +00. Asymptotically, they should obey 



dy f I 'dx(+oo) = dyb/dx(+oo) = \J (q/co) 2 — 1. As for the critical finger, the solution of Eq. 
(PT|) for the front interface can be obtained for any value of the ratio U = q/co > 1 and is 
given by 
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(23) 



v 2 + 1 - 2Uv_ 

with U + VU 2 - 1 < f < +oo. On the contrary, Eq.(p2|) for the back interface can be solved 
with the correct boundary condition only if B is chosen appropriately for each value of U. 
Several obtained shapes are shown in Fig .|. In ref. pi[ , the solution of ( p2|) was computed 
in such a way for several values of U close to one and it was found that c t /c extrapolates 
linearly to one when B — » B c , 

B — B r 



^ = 1 

Co 



K 



(24) 



with the constant K ~ .63. 

We show how the result 
around the critical finger |}| 



can be derived by analyzing perturbatively Eqs.(f21~1, p2|) 
For |ct/co — 1 1 1, the front interface of the retracting finger 
Uf is close to the front interface of the critical finger Yf on distances of the finger tip small 
compared to (q/co — 1) -1 ^ 2 . In this region, we linearize Eq.(pT], |22|) around the critical finger 
shape as yf — Yf + dt/f, Ub = Yb + 5yb- The corrections Syf, 5yb obey the inhomogeneous 
linear equations 
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[B c — - 5B] [Y f (x) - Y b (x)} - B c 8y f 
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(26) 



We have introduced SB = B — B c , 5c t = c t — c and the linear operators Cf, Lb which are 
given by 
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with, 
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(29) 



The boundary values at the tip are 5y/(0) = fo/&(0) = 0. For the derivatives, one obtains 
using the asymptotic behavior Yf(x) ~ — Y(,(x) ~ y2x near x = 0, 
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(30) 
(31) 



As before, Eq. 
of 



can be integrated and one obtains 6yf = rjiSc t /c where 771 is the solution 



'dYj_ 

dx 



(32) 



such that 7/i(0) = 0, 77^ (0) = 1/3. When x — > +00, 771 grows like x 2 /6. The situation 
is different for Eq . (^6|) . For large x, Cb reduces to d 2 /dx 2 — B c . So, in general 5yi> grows 
exponentially like exp(v^B^x) on distances of order one much smaller than the region where 
the linearized equation (p6|) is valid. It is only when <5q/co is related in a particular way to 
5B that the exponential growth is absent and that Syb can grow algebraically like Syj, as it 
should. In order to determine the relation between 5c t /co and SB that should be imposed, 
we find it convenient to introduce the zero mode £{x) of the adjoint C[ which vanishes 
(exponentially) at infinity, 

■72 < 



4(0 = ^ + - = 0, e(+oo) = 



(33) 

where the functions a(x) and b(x) are defined by Eq. (|29|) . ^ is uniquely defined up to a 
global normalization. A local analysis shows that £ automatically vanishes at x = and 
that it tends linearly to zero when x — > 0. For definiteness, we normalize £(x) so that its 
maximum value is equal to 1. A graph of £ is shown on Fig. ^. We now multiply both sides 
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of Eq.(|26"D by £(x) and integrate over x from x\ to xi- Integration by parts gives for the 
l.h.s., 
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dx£(x) C b {5y b ) 



dby b d^ 

~dx ~ dx ~ a ( x >^ x > 5yb 



X2 



X'l 



+ [ X2 dx5y b (x)£t(0 (34) 



The integral on the l.h.s of ( |34"|) vanishes since C\ annihilates £ (|33|) . Moreover, when 
X\ — > and £2 — > +00 the boundary terms also vanish when 5y b satisfies the correct 
boundary condition. Terms at x — +00 vanish when Sy b grows algebraically since £(x) 
vanishes exponentially. There is no contribution at zero since 5y b and £(x) vanish linearly 
which compensate for the singular behavior of a(x) = —3/(2x) + • • •. Therefore, the r.h.s. 



of ( p6[) has to satisfy the solvability condition, 

\h + B c (-I 2 + I 3 )] -5BI 3 = 



Co 



(35) 



where the constants Ji, 13 are given by the following integrals which have been numerically 
evaluated, 
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2.771 



3.814 



7.708 



(36) 



Eq.(p5|) shows that the tangential velocity of the retracting finger tip depends linearly on 
the departure of B from B c as stated in Eq. (p|). The proportionality constant K is in 



excellent agreement with the value obtained by numerical extrapolation in ref. |34 



K = B c + (h-B c I 2 )/I 3 ~ .630 



(37) 



Note that the values (p6|) of the integrals depend on the normalization of £ but that the 
expression of the physical constant K appear as a ratio of such integrals and is thus in- 
dependent of this (arbitrarily chosen) normalization. It is also important to remark that 
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Eq. (^) shows that the tangential tip velocity is an appropriate quantity for a perturbative 
calculation around dR since it has a smooth behavior when dR is crossed. This should 
be contrasted with the retracting tip velocity which decreases like the square root of the 
distance to dR on the retraction side and does not appear to have a simple continuation on 
the spiral side of dR. 

C. Steadily rotating spirals 

For B = e/(a 2 r e A 3 ) < B c , steady spiral waves exist. Their tip rotates around a circular 
core R at a constant tangential tip velocity c t = UiRq. When B — > B c , R diverges, c t — > c 
and the tip of the spiral becomes closer and closer to a critical finger. In this subsection, 
we determine Ro and uu% as a function of B ("the excitability of the medium"). We begin 
by considering spirals of radius large enough so that one can neglect the disturbance of the 
medium due to the spiral previous passage. In this case, the front interface can be assumed 
to propagate in the medium rest state. The spiral shape is analyzed by decomposing it into 
three overlapping regions where different approximations can be performed. Close to the tip, 
on distances of order R tip = e/c , the curvature of the tip trajectory can be neglected and 
a transposition of the analysis of the previous subsection shows that the tangential velocity 
is linearly related to 5B = B — B c by Eq. (^4j) , namely 5c t = 5B/K (both sides being 
negative now). Far from the tip, it is the effect of the interface curvature on the normal 
velocity (Eq. (§)) which can be neglected. The normal velocity can be taken constant, equal 
to Co and the spiral shape is then simply determined. These two approximate descriptions 
match at a distance of order £ from the spiral tip in an intermediate region where the 
interface is almost normal to the tip circle of rotation and the interface curvature is small. 
The intermediate scale I appears as the balance between two effects. On one hand, the tip 
tangential velocity is smaller than cq by about for purely kinematical reasons so that 
| <5c* | ~ cq£/Rq. On the other hand, i is the distance where curvature effects become small 
enough to be comparable to this velocity drop. At a distance £ from its tip, the critical finger 
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curvature is of the order 2R tip /£ 2 . This provides the alternative estimate \Sc t \ ~ eRu p /£ 2 . 
Comparing both expressions and remembering that R tip = e/c gives i ~ (R 2 ip RoY^ 3 and 
\Sct/co\ ~ {Rap/ Rq) 2 ^ . The detailed analysis reported below replaces this simple order of 
magnitude estimate by the precise asymptotic relation, 

S^i^ b (^A 2/3 (38 ) 



Cq V Rq J 



where the numerical constant b is obtained from the first zero a\ of Airy function Ai [42 



b = — 2 1 / 3 ai ~ 2.946 ||36|| . Comparing ( j38l) with ([24]) determines the frequencies and core 
radii of steady spirals near the line A c (e), 



1. Front interface 

We first consider the front interface and assume that the spiral propagates in the medium 
rest state (this is of course justified only if the spiral period is long enough and it requires 
to be sufficiently close to the line OR). The equation for the front spiral interface is then 



identical to the classic equation for the growth of screw dislocations on crystal surfaces [40 
For a steady rotation at frequency u)\ in a counterclockwise direction, Eq.(|) gives for the 
front spiral interface in polar coordinates (r, 6), 



ruji — c Q 1 + r— — + e 



dr 



,Wf + *,1; 2 1 (40) 



V 



dr 



with the boundary condition at infinity d6f/dr — > — oj\/cq. Rescaling coordinates makes 
it clear that Eq. ([4~0"l) depends on the single dimensionless parameter Q = LO\ejc\. For 
< Q < .331, it is found that d6f/dr — > +oo at r = Rq when Eq. (f40|) is integrated 
from r = +oo. Rq is the location of the spiral tip and is found to increase from to +oo 
when Q decreases from .331 to 0. The limit Q — > .331 has been considered in previous works 



3~l] , |2~9"|j . We focus here on the other limit Q — > where the excited region width (~ e/co) 
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becomes much smaller than the core radius (~ ui/cq). In this limit, the front spiral interface 
can be separated into three distinct regions: 

- Far from the spiral tip (the outer region), the interface scale of variation is cq/ui. In- 
troducing the rescaled coordinates r = zcq/uji shows that the terms involving the interface 
curvature are multiplied by the small parameter Q. Neglecting them, Eq. ( [TOD becomes 

d0 ou t 



dz 



Vz^l (41) 



which is of course integrable. This first approximation breaks down near z = 1, where the 
solution of Eq. ( f4T|) has a fast variation on the z-scale and the formally negligible terms are 
important. 

-Close to the spiral tip, Eq. ( |40|) can be simplified in a different way. One can introduce the 
radial distance x from the tip circle of rotation measured in unit of the tip radius such that 
r = R + e/co x and the tangential displacement ye/co = RoOf- At lowest order in e/ (co-Ro), 
Eq. ( f40|) becomes identical in these variables to Eq. ( |i8|) for the front interface of a critical 
finger. 

- These two different descriptions do not directly match. The transition occurs in an in- 
termediate region where the interface curvature is small and the interface tangent almost 
radial. We thus assume (and check a posteriori) that dy/dx is small and expand the square 
root and denominator in Eq. (0). This gives 

c *~ c ° | 6 x - 1 ( d v) 2 | d 'y (42) 
Co co-Ro 2 \dx J dx 2 

where the tangential tip velocity q = Rquj has been introduced. The different terms of 
Eq. (|42| ) are of comparable magnitude for x ~ {RqCq/ t) 1 ^ , dy/dx ~ (e/ ' R$cq) 1 ^ and c t /co — 
1 ~ (e/i?o c o) 2 ^ 3 - in the limit R tip = e/c Ro, this justifies the expansion leading to 
Eq. ( P^) and the neglect of higher order terms. Introducing the rescaled variable £ = 
a;(e/(2c 0j R )) 1/3 , Eq. (g2j) becomes 

ld 2 y 1 / dy\ 2 

2de + i\-di) =e + ai (43) 
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where we have defined 

a 1 = 2- 1 /3^) 2/3 ^ (44) 

The Riccati equation ( |4*3| ) can be transformed into the linear Airy equation. Matching with 
(flT|) imposes that dy/dl; < when £ — > +00. This imposes that the Airy function decreases 
at infinity and be proportional to Ai [f|2|. It gives dy/dl; = 2 + ai)/Ai(£ + ai). Using 
the asymptotic behavior g| Az(f) ~ l/27r" 1/2 £- 1/4 exp(-2/3£ 3/2 ), one indeed checks that 
the obtained large £ behavior dy/dl; ~ — 2^/£ coincides with the behavior of (|4l| ) near z = 1. 
Matching with the tip region requires that the small £ behavior of dy/dl; coincides with the 
asymptotic behavior of ( |I~8"D when x — > +00, namely dy/dx ~ 2/x 2 . This requires a\ to be a 
zero of Ai. Since ?/(£) should be well-defined for all real positive £ it is necessarily the first 
one ai = —2.3381 • • • |4^]. Comparing with the definition (|4]) of a\ directly leads to the 
relation (BSD. 



The relation lui(Rq), numerically determined in j43|, was approximately obtained in 
]1||2"8"| ] by assuming a radial departure of the front interface imposed on a fictitious inner 



radius Ro- This boundary condition is equivalent to requiring that 9 be maximum at Rq . 
It is worth noting that, in the present limit, it would simply amount to replace the exact 
value of the constant a\ by the location of the maximum of Ai namely a[ = —1.01879 • ■ • 
Correspondingly, this would replace the exact value b ~ 2.946 in Eq . (|38f) by b' ~ 1.283. 



2. The back interface 
The equation for the back interface reads in polar coordinates, 

— h^M^T+i**^) (45> 

As for the front interface, we proceed by separately analyzing three regions. We consider 
first the tip region which plays here the dominant role. Introducing as before the coordinates 
x and y such that r = Rq + e/c x, ye/c = RoOb, Eq. ( T45| ) becomes at lowest order in e/c Ro 
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identical to Eq. (|2~2"|) describing the back interface of retracting fingers (except that now 
c t < Co and B < B c . As in this previous case, requiring that the back interface does not 
diverge exponentially from the front interface relates q/co to B. For B close to B c , one can 
linearize around the critical finger and follow the previous analysis ( |25| - |3TD which leads to 
(f2~3|). The comparison of ( |2~3| ) and fl5gp gives the expression ( |55| ) for the spiral core radius 



and frequency of rotation function of B. 

As one moves away from the tip, the back Yj, of the critical finger relaxes exponentially 
toward Yf — 2/B c on the scale of the finger tip width. The equations describing the back 
and front interfaces are thus essentially identical in the intermediate and far regions and the 
analysis of subsection [IV C 1| applies as well to the back equation. 



D. Steady self-interacting spirals 

The analysis of the previous section applies when the spiral period T = 2ir/&i is long 
enough compared to the recovery time constant tr so that the front interface can be assumed 
to propagate in the medium rest state. This applies for A sufficiently close to A c (e) but as 
the medium excitability increases the spiral radius decreases and the front interface begins 
to feel the medium disturbance due to the spiral previous passages. This eventually leads 
to spiral meander as we show in the next section. As a preliminary step, we analyze here 
the influence of this medium modification on the steady spiral parameters Rq and u\. 

The concentration of the controller v on the front and back interfaces follow from Eq. (|7|) 
and (|8|). For a spiral rotating steadily at frequency iO\ in a counterclockwise direction, they 
are given by, 

Vf(r) —v + Svf(r) (46) 
Vb(r) = „ (r) + <W-«k(r) (47) 
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Near the line A c (e), iO\ tends to zero, Svf(r) becomes negligible and the concentration of 
v on the front interface can be taken equal to vq as done in the previous subsection. This 
approximation is justified as long as 8v/ induces a change in the front velocity which is 
negligible compared to the difference c t — Co between the tip velocity and Cq. That is, for 
exp(— 2tiR /c t r ) <C 5c t /co or using the estimate (|55|), exp(—2irR tip /c (bK/B c — B) 3 ^ 2 ) <C 
B c — B. Therefore, one can neglect the perturbation of the medium as long as B c — B <C e 2 / 9 
(up to logarithmic corrections) or equivalently A — A c <C e 5 ^ 9 . The results (^) are modified 
when 5vf(r) becomes comparable to 5c±. The transition regime where 8vf(r) is still small 
can be analyzed along the lines of the previous subsections 
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In the tip region, it is useful to introduce as previously the coordinates x and y, with 
r = R + e/c x,ye/c = R 9. At lowest order in e/(c -Ro) 5 Eq. (|49|) becomes 
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We are interested in the parameter region where aSvf(r)/cQ is of the same order as Sc t /cQ. 
As found above, this happens when the spiral period is large but only logarithmically in e. 
This allows to expand the exponential in fl48|) and to obtain the expression of the medium 
perturbation as a function of the critical finger shape, 



-5v f (r) = B c [Y f (x) - Y b (x)\ exp(-^) 

Co C T R 



(52) 



We expand the spiral front in the tip region around the critical finger shape as yf(x) 
Yf(x) + 5yf(x). The correction Syf(x) obeys the equation 
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(53) 



where the linear operator Cj is defined by (|27|) . Syf can be expressed as 



27 



5c t 



5y f = — 771 + B c e r) v 
c 



(54) 



where r/i is defined in ([32|). 77^0 obeys 

C f ( Vvfi ) = (Y f (x)-Y b (x)) 
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(55) 



with the boundary conditions ?7„,o(0) = 0, ^p(O) = 2/3. 

In the same way, we obtain, in the tip region, the lowest order correction 5y b to the back 
interface of the critical finger Y b (x) 
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(56) 



where the linear operator C b is defined by Eq. (^). Eq. (^) is similar to Eq. (|26|) and 
can be analyzed in the same way. The solvability condition that should be obeyed in order 
for 5yb not to diverge exponentially as x —>■ +00 is found by integrating both members of 
( j56|) with the zero mode £(x) of the adjoint of C b . This gives the following generalization of 
Eq. (H 



K — = 025 + B c J exp( 
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where K ~ .630 (Eq. (^))and the constant J is 

J = 1 + B C ^ ~= 1.872 
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-Z3 is defined by ( p6|) and I v $ is given in terms of £(x) (|33D and 77^0 
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To complete the analysis, it remains to match the tip region to the outer part of the spiral. 
As one moves away from the tip, the finger width Yf(x) — Y b (x) relaxes exponentially toward 
its asymptotic value 2/B c on the scale of the tip region. Therefore in the intermediate and 
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far region, a8vf(r)/co fl52"D is equal to lowest order to 2exp(— 2-kRq/ cqTr) and the matching 
equation becomes instead of (pE2] ) 

ct - cp 2vri?o e I ( dyV d 2 y 

+ 2 exp + ——x = - — + — 60 

Co Cor^ co-rto 2 \ax J ax z 

Matching with the tip region gives in a similar way 

'c R \ 2/3 



2 -l/3 



c t - c , 2ixRq , 

+ 2 exp( 



a x ~ -2.338- •• (61) 
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The difference between Eq. ([Jl]) and ([38|) is simply that c t is not compared to the velocity 
of a single planar pulse but to the velocity of a train of pulses of wavelength 2nRo (the 
asymptotic wavelength of the spiral to lowest order). Comparing (|57|) and (^1|) determines 
the radius R$ of a steadily rotating spiral as a function of the medium characteristics B, 



B c -B = Kb[ — — + (S c J + 2it) exp( -) (62) 

The medium disturbance due to the spiral previous passage has two distinct effects which 
are comparable to lowest order: 

-the medium "excitability" is reduced in the tip region which modify the tip velocity 
(Eq. ©), 

-the tip velocity should be compared to the velocity of a periodic train of planar waves which 
is slower than the velocity of a single planar pulse. 



V. DERIVATION OF KINEMATIC THEORY 

We consider now the spiral dynamics in the vicinity of the line A c (e) (for e<l). In this 
limit, several simplifying features made the previous analysis of the steady states possible. 
These still hold when one is interested in an unsteady motion taking place on a time scale 
comparable to the steady rotation period which is long compared to the time scales of the 
internal modes of the wave tip: 

- i) the dominant effect which shapes both the steady spirals and retracting fingers tips is 
the curvature dependence of the normal velocity. As a consequence, the shape of a wave tip 
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is close to a critical finger up to a distance I from the tip where the curvature effects have 
become small enough to be comparable to the velocity difference between the tip and planar 
front velocity, namely when Co — ct ~ eRtip/£ 2 (where we have evaluated the curvature 
—d 2 Yf/dx 2 ~ Rtip/V 2 at x ~ t using the asymptotic behavior Yf/R tip ~ \n(x/Ru p ) for 



x/Rup ^> 1). This yields the relation t ~ R tip /^1 — c t /co that remains also true in the 
unsteady case. The motion of this 'solid' shape can be determined from the knowledge of 
its instantaneous tangential velocity q and of its instantaneous rotation rate u obtained by 
extending our previous analysis of the steady states, 

-ii) the tangential velocity q depends on the 'average' concentration of the controller v in 
the vicinity of the tip. The precise definition of the average is obtained by using a solvability 
condition which generalizes Eq.(|35|) and fl57|). 

-iii) a tangential velocity c t smaller than the asymptotic normal velocity Cq of the wave gives 
rise to a rotation of the solid tip at a rate u which can be estimated as in the steady case. 
As said above, the wave tip has a solid character (i.e. is close to a critical finger shape) up 



to a distance £ ~ Ru v / y 1 — q/co from the tip. Since kinematics requires that u£ ~ Co — c t , 
one obtains for the rotation rate lu ~ cq/Ru p {1 — co/c t ) 3 ^ 2 . As shown in Appendix (B), this 
relation is established on the time scale ~ Rt ip /(c — c t ) much shorter than the steady 
rotation period. Therefore, on this latter slow time scale, the slowly varying rotation rate is 
linked in an adiabatic manner to the slowly varying tangential velocity by the same relation 



Eq. (|38D or (pTj) which relates the steady state frequency to the tip velocity. 

We begin our analysis by considering the kinematics of the wave tip motion. We then 
compute the tangential velocity of the tip as a function of the concentration of the controller 
v in the medium left by previous passages of the wave. As a result, we obtain an ordinary 
differential equation with delay which describes the motion of the wave tip. An analysis of 
this equation at the linear and the weakly non-linear levels determines the characteristics of 
the meandering instability near threshold in the weak excitability limit. 
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A. Parameterization of the wave tip motion 



We use polar coordinates (r, 9) with the origin at the center of the circular steady spiral 
core. The wave tip motion is determined by its tangential velocity c t ({v}), a functional 
of the (space and time dependent) controller concentration which will be computed in the 
next subsection, and by the rotation rate of the shape (or, equivalently by the radius of 
curvature of the tip trajectory). We use a complex notation z{t) = r(t) exp(z#) to denote 
the tip position. Then, the tip velocity is \z\ and the shape rotation rate Im(zz/\z\ 2 ) where 
time differentiation is denoted by a dot and z is the complex conjugate of z. The tip motion 
is thus determined by the two equations 

\z\ = ct({v}) (63) 
Im(zz/\z\ 2 ) = c t ({v})/R t [ Ct ({v})} (64) 

where at this stage q({w}) can be thought of as a given function of time. The instantaneous 
radius of rotation Ri is a function of c t ({v}) given to lowest order in the interaction parameter 
by Eq. (0), 

'c 0j RA 2/3 



ct - c , 2-nRo , 

+ 2 exp( 



(65) 



Co C T R 

We will actually find that the meander threshold occurs before a significant modification of 
the steady state radius by the interaction so that Eq. fl65|) can be replaced by the simpler 
Eq. © 

RiM^-i-^-V* (66) 
c o \ Co — c t ) 

We consider the motion of a spiral tip which is displaced from its steady state position 
z = (Ro + eq(t) / co)e tullt+ ^ (see Fig. |3|). We restrict ourselves to displacements of the tip 
which are comparable to the tip radius of curvature e/co ( i.e. q(t) ~ 1) and therefore small 
compared to the radius of the steady core i?o- As a consequence of the tip displacement, 
the controller concentration and thus Ct({v }) and Ri depart slightly from their steady-state 
values, c t ({v }) = c° t + Sc g ({v }), Ri = R + 5Ri. 
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We assume (and will check a posteriori) that the time scale of the unsteady motion is of 
the order of the steady state period T = 2ir/u\. We expand Eq. flB5| ) and ( |6lD in the small 
parameter e/(co-Ro) and keep only the dominant terms, \z\ = luiRq(1 + ip/ui + ge/(co-Ro) + 
■ ■ •), Im(zz /\z\ 2 ) = Ui(l + ip/uji — eq/ {cqUj\Rq) + ■••). Eq. (|63|) gives therefore at lowest 
order in e/(c i?o) 



■U) X q 



Rq CoRo 

which shows that ip/ui ~ e/(c -Ro)- Using this scaling, Eq. 

e q 



(67) 



becomes at lowest order 



_ c ° ^ i 

R 



* Rl 



(68) 



We obtain the equation for the radial motion of the tip by substituting in (|68"D the expression 



q + ufq = Luj5RiC /e 



(69) 



Finally, it is convenient to use the tip angular position 9 = io\t + ip(t) instead of time. To 
lowest order in q/R, this simply gives 

cvR'Acf 



d?_q 

de 



2 +q 



(70) 



with 



cpijM = _3_ r cpRp 
e 2b Co V e 



5/3 



(71) 



from a differentiation of Eq. (p5|) . In order to have a closed equation for g, it remains to 
express Sc q in terms of the previous positions of the wave. We now proceed to this task. 

B. Computation of the tangential tip velocity for self-interacting spirals 

We consider successive passages of the wave tip by the angular position 6. The successive 
radial displacements of the tip are • • • , eq{6 — 2n)/c , eq(d)/co , eq(9 + 2k)/cq , • • •. Let us 
consider the passage at the position Rq + eq(9)/co in the Cartesian coordinate system (x, y) 
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attached to the wave tip (see Fig. |^) in which we choose to measure lengths in unit of e/co . 
The controller concentration 9) in the medium just ahead of the front interface is related 
by the controller recovery kinetics Eq.(^) to the controller concentration left just behind the 
back interface Vb(x;9 — 2tt) at the previous passage. At dominant order in e/(coRo), one 
can neglect the tip width compared to the core perimeter and the time interval between two 
passages of the spiral by the same angular position can be taken equal to the steady spiral 
period T , 

v f {x- 6)-v = exp(-^) [v b (x + q{9) - q{9 - 2tt); 9 - 2tt) - v ] (72) 

In (j72|), note that Vf(x;9) is related to Vb(x + q(9) — q(9 — 2n); 9 — 2n) since the argument in 
Vf refers to a frame attached to the wave tip with origin at Rq + eq(9)/co whereas the origin 
of the coordinate for Vb is at R + eq(9 — 27r)/c . 

The controller concentration Vf and Vb at the same passage are also simply related by 
the controller production equation in the excited region (|7|), 

v b (x; 9) — Vf(x; 9) + e— (73) 

C QT e 

Iterating back in time ( [72] ) and ([T3|), we see that f/(x; 9) depends in principle on the positions 
of the tip, at all previous passages by the angular position 9. However, the memory of the 
position q(9 — n2n) is suppressed by the n-th power of the small parameter exp(— T /tr). 
Therefore, to dominant order the controller concentration only depends on the position of 
the tip at the previous passage 

„,(*; 9) = v + e e~^ R vj? + <M Z + d(6)) ^ + d{$)) (?4) 

where we have defined the relative displacement of the tip between its two passages d{9) = 
q(9) — q{9 — 2tt). is the usual Heavyside step function, B(x) = for x < and Q(x) = 1 
otherwise. Eq. ( |74|) determines the controller concentration on the front interface at 9 as a 
function of q{9) — q{9 — 2tt). It is now an easy task to generalize the previous computations 
and obtain the tip tangential velocity corresponding to this concentration. 

33 



As for steady interacting spirals, we obtain for the front interface in the tip region, 



d 2 y/ _ ct_ 
dx 2 Co 



1 + 



dx 



a: 



1 (v f (x;6) -v ) 

c 



l + { djH~ 



dx 



3/2 



(75) 



The only difference with Eq. ([51]) is that Svj = Vf(x;9) — vq is now given by Eq. ( [T4| 
Expanding Eq. ( [75] ) around the critical finger shape, yf = Yf + 5yf, one obtains as before 

2 



Co 



dK 



/ 



dx 



+ B c e- To/TR [Y f (x + d{9)) - Y b {x + d{9))\ 9(x + d{9)) 



'dYj_ 
dx 



3/2 



(76) 



with the solution 



Sc 

<% = — % + B c exp(-T / t r ) r] vAe) 
Co 



(77) 



The linear operator Cf is defined by (|27|), 771 is defined in (|32| ) and is the solution of 

3/2 



= [Yf{x + d) - F b (x + d)] 9(x + d) 



1 + 



dF, 



/ 



dx 



(78) 



which generalizes Eq. (|55|), with the boundary conditions at x = 0, ^,d(0) = 0, rj Vt d{x) 



x/2[Y f (d) - Y b {d)} 6(d) for x < 1. 
Similarly, the back interface equation in the tip region is 



d 2 y& _ Ct_ 
dx 2 c 



1 + (^) 2 
dx 



1 MX; 9) - v ) —(yf(x) - y b {x)) 

Co C t C Q T e 



1 + (^) 2 

dx 



3/2 



(79) 



After linearization around the back interface of the critical finger, y&(x) = Yf,(x) + Syb(x), 
one obtains for the correction 



Set 
c 



1 + 



'dny 

dx J 



fir 

[B c — - SB][Y f {x) - Y b {x)} - B c e- ToTR [Y f {x + d{9)) - Y b {x + d(6))]G(x + d{9)) 



c 



--B c Sy f 



1 + 



dn 

dx 



3/2 



(80) 
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Multiplying both sides of (p0|) by the zero-mode of the adjoint of C b and integrating 
from x = to +00, gives 

Set 



Co 



[h + S c (/ 3 - J 2 )] = 6B h + B c e- Ta ' TR [h m + BJ V 



d(9)\ 



(81) 



where the definite integrals h, h have been defined in ( pop and I^d, I V: d are given by 

3/2 

r+x / <IY,.\- 



dx^(x) [Y f (x + d)- Y b (x + d)] G(x + d) 

3/2 



1 



dx 



dx£{x)r] V)d 



1 + 



d\l 

dx 



(82) 



Finally, this gives the tangential tip velocity as a function of the tip displacement 



5c t 5B B c 



-T /t r 



[J + F(q(9)-q(9-2n))} 



(83) 



c K K 

where the constants K ~ .630 and J ~ 1.872 are defined in Eq. ([37]) and Eq. (|58|). The 



function = [(-^3^ + B c I v>d )/I 3 — J] vanishes at d = and is plotted in Fig. || [44 . 

Comparing Eq.(|83|) with Eq.(^) for the steady case shows that the change in tangential 
velocity due to the tip displacement is 



5, 



Sc q = c -^ e- To/TR F{q{6) - q{6 - 2tt)) 



(84) 



C. Computation of the tangential velocity in other cases 

We conclude this section by emphasizing that, although the present kinematic theory is 
quite general, the precise expression for the tangential tip velocity that is to be used in in 
conjunction with Eq. [70] depends on the application at hand. For example, Eq. ^4| above 
is valid for self-interacting spirals without external forcing and is therefore perfectly suited 
to analyze meander in the next section, or interacting multi-armed spirals with a minor 
modification given in section | V 1 1 1| . For the non- self-interacting spiral with an excitability 
that varies slowly in space or time (section VILA), one can use directly the results for steady- 
state rotation (Eq. P4|) , whereas under the action of an external field (section VII. B) one 
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needs to compute a different expression for the tangential velocity. The general procedure, 
however, is clear. In each case, the tangential velocity depends on the average controller 
variable in the tip region and can be computed from a solvability condition. 



VI. MEANDER 

In this section we analyze the classic meandering instability in a linear and nonlinear 
regime. Substitution of (jS4]) in ([70|) expresses the r.h.s. of (|70|) as a function of q{6)—q(6—2-n) 
and provides the differential equation with delay governing the tip motion 



dO 



-jj± + q = mF(q(6)-q(9-2n)). (85) 



The parameter m is given by 

m = H (^) 6 " exp(-T /T R ). (86) 

Values of m of order unity are reached when (R co/e) 5 ^ 3 exp(— T /r/j) ~ 0(1). In 
this parameter regime, one can use the simple formula (|39| ) to estimate the spiral pa- 
rameters since in ( |62|) the correction term (the second term on the r.h.s) is of order 
(i?o c o/ e ) 2//3 exp(— 2ttR /c t r ) ~ 0(1) compared to the first term on the r.h.s. and therefore 
smaller by e/ (RqCq). This provides the explicit expression of m in terms of the parameter B 
which characterizes the medium 



3B c (bKf/ 2 
m = 2{B C - B)W 6XP 



2vre / bK \ 3/2 ' 



c 2 r R \B C -B / 



(87) 



A. Linear stability analysis and instability criterion 

We begin by studying the linear stability of Eq. ( |85| ) around q=0 that is, the linear 
stability of steady rotation. For q 1, one obtains 

^ + q = a (q(6)-q(6- 2*)). (88) 
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where we have introduced a = mF'(0). Seeking q under the form q = Aexp(a9) gives the 
eigenvalue equation 

a 2 + 1 = a [1 - exp(-27nr)] (89) 

For any a, a± = ±i are isolated solutions of (|89|). They simply correspond to the two 
translation modes of the spiral : for a steady spiral which is slightly displaced from the 
origin and centered at (x ,y ) with x <C R, yo <C R, the distance of the wavetip to the 
origin varies sinusoidally as q = \zq + Ro exp(i9)\ — Ro = xq cos 9 + yosin 9. 

The other solutions of (^) vary with a. For small a > 0, the r.h.s. of ( p9|) is comparable 
to its l.h.s. only if the real part of o is large and negative, that is Re[o) ~ — l/2ir ln(a). 
Therefore, for small a, all eigenvalues (different from the two translation modes) have a 
negative real part and the steady rotation is stable. As a is increased, the eigenvalues moves 
continuously in the complex plane. An instability occurs when the real part of some of 
them traverses zero and becomes positive. This happens at the critical value a = a c where 
Eq. ( |89D has a purely imaginary root a = ifl, namely for 

a c [1 - cos(27rfi)] = 1 - tt 2 (90) 
a c sm{2nQ) = (91) 

Eq. (|9~lD requires that Q be a half integer. Eq. (|90|) can therefore be rewritten as 1 — Q 2 = 
a c [l — (— l) 2n ], the only solution of which is, for a c > 0, = ±1/2, a c = 3/8. 

We therefore conclude that for < a < a c all eigenvalues different from a± have a 
negative real part. As a increases past a c = 3/8, a couple of eigenvalues traverses the 
imaginary axis and acquires a positive real part. The value a = a c is thus the threshold of a 
Hopf bifurcation and corresponds to the meander onset with a frequency ratio at threshold 
u)iju)\ = 1/2. This ratio is consistent with the extrapolation to infinite core radius of 



numerical simulation results as shown in Fig. 10 



It is interesting to note that as a is further increased, the frequency of the two linearly 
unstable modes decreases and the two unstable eigenvalues become purely real for a > a r 
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(a r is simply determined as the value of a for which Eq. ( |89"D has a doubly degenerate 
root, «fl = a r /n exp(27rov) with af + 1 — a/ir (exp(27ra r — 1) which gives a r — .375 and 
a r ~ 1.260). This may explain why a previous analysis performed at small e, but away 
from dM ||29|| , yielded only real unstable modes instead of complex conjugate eigenvalues as 
expected from a Hopf bifurcation. 



Given the expression (|87|) of the constant m, the criterion for meander onset a c = 
m -F'(O) = 1/2 implies that, for small e, the meander boundary A m (e) lies close in the (e, A) 
plane to the critical finger boundary A c (e) (see Fig. I) with A c -A m ~ e 5/9 / ln 2/3 (e 5/9 /T'(°))- 

In the pure sharp boundary description with no diffusion of the controller v field, the 
behavior of the function F is nonanalytic at short distance F(d) ~ —.576c? ln(|d|) for d <C 1 
as shown in Appendix 0. This implies that, in this context, the onset of meander occurs right 
at the critical finger boundary. However, when one starts from the full reaction-diffusion 
equation (|],[D, the interface has a finite width of the order of e. This eliminates any short 
distance nonanalyticity and cut-off the divergence of F'{d) at d ~ A ~ e 1 / 3 . This gives the 
estimate F'(0) ln(e) and A c - A m ~ e 5 / 9 /| ln(e)| 2 / 3 . 

The non-analyticity of F also disappears if the slow field v diffuses that is, if instead of 
([]) one has 

^ = 7£V 2 H5(fJ ± (v),t)) in V ± (92) 

For a sufficiently small diffusion constant 7, one can neglect entirely the diffusion in the ex- 
cited region T> + and consider only a radial diffusion of v in T>~ . The controller concentration 
on the front spiral interface is then a smoothened version of 



v f -v = / e -W)-*?l* D {Yf (x') - Y b (x')) (93) 

v b = v f + e (Y f (x) - Y b (x)) /(c c t r E ) (94) 



The finite diffusion length Id = ^/AeyTo co/e removes the short distances analyticity and 
gives a finite first derivative to F at the origin which decreases with increasing ^as plotted 
in Fig. H (see also Appendix 0). This decrease of stability with a decrease of £d qualitatively 
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agrees with the numerical results of ||45|| . Of course, diffusion controls the stability only if 
Id is much larger than the interface width (or the width of the tip boundary layer). When 
it is much smaller, stability is controlled by finite interface width effects as discussed above. 



When the two effects have comparable magnitude, the numerical results of [f5J suggest that 
more complex stability diagrams are possible (i.e. there is a region of reentrant stability). 
It would be interesting to see if this could be explained by a more complete computation of 
F taking into account both finite interface effect and diffusion of v. 

We conclude this subsection with a simple interpretation of the obtained results. The 
existence and magnitude of the instability threshold can be understood by considering a 
displacement of the wave tip by a small distance d towards the outside of its steady circular 
trajectory. Since the outside of the core is slightly less excitable than the inside, this outward 
displacement will cause the spiral tip to propagate in a less excitable medium and to rotate 
on a new larger radius Ri = Rq + SR^ > Rq. The fact that 5Ri > by itself is not sufficient to 
create an instability. It is only if 5Ri is larger than ~ d that the displacement of the wave tip 
can be amplified and meander can appear. The excitability change due to the displacement d 
is \SB\ ~ dj R tip e~ T °/ TR where the exponential factor simply reflects the global attenuation of 
excitability variations between two passages of the wave. This excitability change leads to a 
variation of the rotation radius SRi ~ dRo/dB SB. Thus, SRi/d ~ (dRo/ dB) / Ru p e~ T °^ TR ~ 
m and the onset of meander occurs for m of order unity in agreement with the above stability 
analysis. The period doubling like character of the unstable motion (i.e. uj 2 /uji = 1/2) can 
also be attributed to the radial gradient of excitability at the edge of the spiral core. A 
wave tip displaced outward from the center at a given passage, will propagate, at its next 
passage, in a medium more excitable than the one produced by steady rotation. This will 
cause the tip to execute this second turn on a smaller radius and thus, to propagate again 
in a less excitable medium and with a larger radius at the next cycle, leading to the period 
doubling behavior. As we shall discuss in section IX, this picture is modified by finite core 
effects that roughly make trajectories of larger radius take a longer time to complete one 
rotation. This effect causes the spiral tip to return sooner inside the core and, in turn, leads 
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iOiju)\ to increase away from 1/2 with decreasing Rq. 



B. Nonlinear dynamics 

We now carry out a standard weakly nonlinear analysis of the wavetip equation of motion 
(jB5|) and show that the bifurcation to meander is supercritical in agreement with existing 
numerical studies of reaction-diffusion models |14| , |15[ . This analysis also allows us to char- 
acterize more precisely the epi-cycle-like trajectories of the wave tip in the large core limit. 
Next, we integrate Eq. ^ numerically and explore the nonlinear regime further away from 
the bifurcation point. 



1. Weakly nonlinear analysis 

To carry out the weakly nonlinear analysis, we first expand the function F on the r.h.s. 
of Eq. B5| up to cubic terms, which yields the equation 



+ q = aAq + T(Aq) 2 - (3(Aq) 3 (95) 
where we have defined 

Aq = q(6) - q(6 - 2tt) (96) 

and the constants 

a = mf'(O) (97) 

T = mF"(0)/2 (98) 

= -mF'"(0)/6 (99) 

In writing (p5|), we have supposed that the non-analyticity of F in the pure sharp boundary 
limit has been taken care of either by taking into account finite interface width effects or by a 
small diffusion of v (e.g. for £ D = 1 one has F'(0) ~ 1.12, F"(0) = 2.81CT 2 , F w (0) ~ -1.1). 
Note however that Eq. (|85|) is well defined even for the non-analytic sharp boundary F . We 
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shall comment in the next subsection on the small amplitude behavior in this case. Eq. 
is valid in a regime where the parameter 



/i = a — a c 



(100) 



which defines the distance above the onset (a c = 3/8) of the meandering instability is small. 
Next, we seek perturbatively for time periodic solutions of Eq. |95| of the form 



Q(0) = % + E A nZ m ™ + CC, 



(101) 



n=l 



where as before Q = uj 2 /uj\ is the ratio of the Hopf frequency at the meander bifurcation 
and the primary angular rotation frequency. Substituting Eq. |101| into Eq. ^ and focusing 
on the first two modes (n = 1 and n = 2), we obtain at once that 



'-n 2 + l)A 1 = a(l - C)A 1 + 2TA 2 A 1 (1 - C)(l - C) 

- 3/3(1 -C) 2 (l-C)^i|^i| 2 
-4fi 2 + l)A 2 = r(l - C) 2 A\ + a(l - t 2 )A 2 

- 2/3(1 -o(i- o(i -mi^i 2 



(102) 



(103) 



where we have defined 



C = e 



(104) 



and (, A n , denote the complex conjugates of (, A n , respectively. Eliminating A 2 between 
the above two relations, and neglecting the terms proportional to (1 — () 2 (1 — O^bl^UI 2 on 
the r.h.s. of Eq. |102| (which can be checked to be of higher order at the end), we obtain that 

, 2 (i-C))i-C) 2 (i-C 2 ) 



n 2 -l + a(l-() + 



2P 



3/3(1 -cr(i-c) 



|Ai| 2 = (105) 



1 - 40 2 - a{\ - C 2 ) 

The condition that the real and imaginary parts of the l.h.s. of the above equation must 
vanish independently provide two independent relations that determine Q and A\. Next, 
expanding Eq. |105| to first order in the frequency shift Q — 1/2, we obtain 

16ir 2 vr 



- 3/4 + 2a + (1 - i2na c )(n - 1/2) - 



24/3 + 



1 + ia c n 







(106) 
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The conditions that the real and imaginary parts of the above equation must vanish lead 
after simple algebraic manipulations to the relations 

\Al\ = ^C~JI (107) 

n - 1/2 = -C2IM (108) 

where C\ and c 2 are constants defined by 

l/ci = 12/3+— , ,\ (109) 

1 H 3 1 + (3tt/8) 2 V ; 

64 Cl r 2 

C2 " 3 [1 + (3vr/8) 2 ] (110) 
Eq. p.03| implies that at leading order in //, 

A 2 = — (111) 

(1/2 — n) (1 + i 3tt/8) V ; 

or, using Eqs. |107| and pL08| , 



\M = Trrp + (3vr/8) 2 (112) 



64r 

In addition, substituting Eq. |101| into Eq. |95|, one obtains for n = that g — 8T\A 



|2 



8rci/i. It is simple to work out that higher order terms in the present expansion must scale 
as A n ~ /i n//2 for n odd and A n ~ yU n ~ 2 for n even. Note that the expansion of (|85|) leading 
to (^5| ) remains justified because Aq(6) vanishes as [i — > even though A2 remains of order 
unity (i.e. A 2 [exp(i2Q6) — exp(i2Q(6 — 2n)] ~ c 2 fi in this limit). 

Let us now examine the meander trajectory of the wave tip. For this purpose it is 
convenient to define the dimensionless coordinate Z = X + iY = Re /Ru p , which is scaled 
by the tip radius Ru p = e/co, and is given by 

Z = X + iY = (p + q)e iie+ ^ (113) 
diP/d9 = -(q-q )/po (114) 

where we have defined the scaled steady-state radius po = RoRtip- We have subtracted the 
^-independent part of q(9) which gives a shift of uj\ of O(go/po) (Eq. (p7|). Since 1/po <C 1, 
we can expand the above relations to first order in if), which yields 
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X + iY = (^p + q-i J(q-q ) dd^j e i{1 ^ o/pQ)e+i ^ (115) 

Since the phase factor i[) corresponds to a translation of the center of rotation, we can set 
■00 — 0, which yields the relation 



X + tY = p e w + J2 



n=l 



(116) 



where the amplitudes A n dictate the meandering motion of the tip. 

Note that in deriving Eq. |116| we have only assumed that q/po is small, such that this 



equation is not restricted to the asymptotic large core limit where Q = 1/2 at the bifurcation. 
In fact, in the weakly excitable limit that is typically accessible in simulation, Q is larger 
than 1/2 at the bifurcation point due to finite core radius corrections ~ 1/po that modify 
Eq. |85] as discussed in section |IX| (see e.g. Eq. ( |165|) ). In this case, the bifurcation is not 
resonant (i.e. 2u 2 ^ uj±), and A 2 ~ p near onset. Eq. |116| implies that in this generic 
case, relevant for usual simulations and experiments, the motion of the tip can be described 
by keeping only the terms proportional to A\ and A\ in Eq. ( |116| ) that is a three-radius 
epi-cycle (or epi-epi-cycle) 

X + iY = p e id + p x e ^-^-^ - P2 e <(i+n)«-H«i > i / 2 ) (117) 

where 0\ is an arbitrary phase, pi ~ JJi, and 

P2 / Pl = (i-n)/(i + n) (us) 

The fact that p 2 /pi vanishes as Q —>■ 1 may provide an explanation for why the meander 
trajectories in simulations of react ion- diffusion models of excitable media have been tradi- 
tionally well fitted by a simple epi-cycle (Eq. |117| with p 2 = 0). In Ref fl5| , it was argued 



that meander trajectories should generally be epi-epi-cycles close to the onset of instabil- 
ity. It was left unexplained, however, why the ratio p 2 /pi turns out to be very small. For 
the simulation of the FN of Ref. JT5] , Q ~ 0.782, in which case Eq. |118| predicts that 
P2/P1 ~ 0.12. This ratio is roughly consistent with the ratio of the amplitudes of the peaks 
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of 1 +Q and 1 — Q in the power spectrum of X(t) in Fig. 4 of [[nj. Here, Fig. illustrates 



that two-radius and three-radius epicycle trajectories are very close even when Q departs 
significantly from unity. Such a small difference is probably hard to resolve experimentally. 
Let us now examine the meander trajectory predicted by Eq. |116| in the asymptotic limit 



where Q = 1/2, which is more difficult to reach in simulation and experiment. The main 
difference in this case is that A2 is 0(1) because the bifurcation is resonant, i.e. A\e e act 
as a periodic drive of the wave tip at the primary frequency Q = 1. Inserting the results of 
the weakly nonlinear analysis, Eqs. |107| - |112| , into Eq. |116| , we obtain that 



X + tY = p o e i0 + Pl e^ 2 -^ - p 2 e m / 2+i0 i + ^^(-^i+tan-^/s) (0 = 1/2) (119) 
where we have defined Q* = 1 — 2Q = 2c2/i, 6± is an arbitrary phase, and 



Pl = 3V^Z, p 2 = Pi/3, p 3 = 3^1 + (3tt/8) 2 /(32 T) (120) 

Consequently, the effect of the resonance when O = 1/2 is to add a slow component of 
motion with frequency Q* ~ p around a circle of radius p% of 0(1). Steady-state rotation 
is approached smoothly when fi — > 0, even though p 3 remains finite, because Q* vanishes 
in this limit. Finally, we note that p 3 diverges as 1/r in the limit r —>■ 0. The tangential 
velocity of the tip around the circle of radius p 3 , however, scales as £l*p3 ~ T and vanishes 
in this limit, which is therefore well-behaved. 

2. Numerical integration of the wavetip equation 



Eq. 85 was integrated numerically using the algorithm described in Appendix A. We 



used both the function F plotted in FigJ^, and the simple analytical form 

F(x) = tanh(x - a) + tanh(a), (121) 

This form has qualitatively the same shape as the calculated function F, which is plotted 
for different 1^ in Fig. ^, and yields a qualitatively similar nonlinear behavior. For this 
reason, all the results presented here are for this simplified form of F defined by Eq. |121 
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for the choice of parameter a = 0.2. As noted earlier, calculated function F is non analytic 
at the origin in the singly diffusive sharp boundary model and behaves as — .57Qq m(|g|). 
When this is used in Eq. (|5|) for the tip motion, as noted previously, a steady rotation is 
unstable for all m (Eq. (p7|)) however small since the slope of F at the origin diverges. The 
growth of the modulation as one moves away from threshold is however much slower than 
in the analytic case, the amplitude of the modulation being of order exp(— cst/m). It is 
interesting to note that requiring this amplitude to be larger than the interface width e, as a 
criterion for meander threshold in a real small-e model, gives m ~ cst/\ ln(e)| quite similarly 
to what was obtained previously by cutting off the slope of F at the scale of the interface 
width. Away from onset, however, this non-analyticity does not modify much the nonlinear 
behavior. For this reason, we shall not treat this case separately. 



The results of the numerical integration of Eq. |85] are illustrated in Figs. |TT| and |12[ We 
have found it convenient to plot q{6) — q{0 — 2tt), instead of q{6) because the latter quantity 
contains a component ~ e that only yields a translation of the center of rotation. We have 
checked that the amplitude of oscillation and the frequency shift of Q from 1/2 increase 
quantitatively for small /i as predicted by the weakly nonlinear analysis. Fig. IT2| shows that 



the oscillations become more nonlinear with increasing distance from the bifurcation point, 
but remain periodic with a frequency close to 1/2. The fact that the frequency is rather 
insensitive to m can be understood by remarking that F (calculated or approximated by 
Eq. |121| with a small) is close to being an odd function of its argument. For F exactly odd 



(r = mF"(0) = 0), the weakly nonlinear analysis of the last section predicts that A n = 
for all n even and that there is no nonlinear frequency shift, i.e. Q = 1/2 for any value of 
fi > 0. One would therefore naturally expect to to find that Q remains close to 1/2, even 
far from onset, when F deviates slightly from an odd function. 

Finally, it is worth noting that hyper-meander (i.e. chaotic meander) is not contained 
in the large core limit. This is consistent with the fact hypermeander has been observed 
numerically in the opposite parameter range of high excitability [|12|| . In this range, the 
shape of the spiral boundary is not constant in time on the scale of Ru p - It therefore seems 
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likely that the dynamics on this scale plays an important role in hyper-meander. 



VII. SPIRAL MOTION UNDER EXTERNAL ACTION 

Motion of spiral waves can be induced by modulating the medium excitability in space 
or time or by adding an external field. It is not difficult to extend the approach of section 
[V|to describe these effects simply and quantitatively in the large core limit. 

A. Variation of the medium excitability 

We consider first the effect of spatial and temporal modulations of the excitability (ob- 
tained by changing A and/or e into space and/or time). Such a modulation will generally 
produce a variation of both the planar front velocity Cq and a variation 5B(z,t) of the pa- 
rameter B characterizing the medium. We assume that this variation is small enough to be 
treated as a perturbation, that 5B(z,t) varies slowly in time (i.e. on the scale of the spiral 
rotation period) and in space (i.e. on the scale of the spiral core) and that B is close enough 
to B c (i.e. dR) so that the spiral self-interaction can be neglected. The radius of curvature 
of the tip trajectory will then depart from its unperturbed value R , Ri = R + 5Ri with 

« = (122) 

Ro 2 B c — B 

and the variation of Cq gives a subdominant contribution for B close to B c . Substituting the 
above expression into Eq. (|69|), we obtain at once 

2 3c R ujl5B(z,t) 

9 + ^=2— -aTTB- (123) 



Integration of Eq. ( |123| ) gives the spiral tip motion resulting from a given space time 



variation of excitability. As a simple illustration, we show that a global periodic variation of 



excitability at the spiral frequency induces a spiral drift |]I5f . When 5B = A cos(u>it+(j)), the 
r.h.s. of Eq. ( [f23| ) is resonant with the natural oscillation modes of the l.h.s., the translation 
modes, and induces their growth 
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g (t) = ^£^AsinM + 0) (124) 

A simple way to understand the motion described by ( |124j ) is to remember that for a 
steady spiral centered close to the origin (compared to the radius of its core), at z = xo+iyo, 
the distance of the wave tip to the origin varies periodically as 

\Rq exp(iuit) + z \ ~ R + Xq cos{u\t) + y sin(uit) (125) 

Comparing the two expressions shows that Eq. ( |124| ) describes a linear drift of the spiral 

^ = 3 c q Rq ^ _ B uit (-i exp(-i0)) (126) 

The drift direction depends on the relative phase between the spiral rotation and the periodic 
modulation of excitability : the spiral drifts perpendicularly to the direction (exp(— icf))) of 
the spiral tip at the maximum excitability viewed from the spiral center. One can note 
that our derivation of ( |126[ ) is simple, but, of course, it breaks down when the spiral center 
is no longer close to the origin and the linearization giving fl6"9|) and thus ( |123| ) becomes 
illegitimate. The remedy is standard: a nicer looking derivation is obtained by introducing 
a the very start of the derivation of ( |69|) the spiral center zq and parameterizing the wave 



tip as z = z + (R + eq(t) / c )e luJlt+ ^^ . The slow variation of Zq with time is obtained by 
requiring that it cancels the secular term on the r.h.s of Eq. (|123|) . 

A time independent excitability which varies slowly in space is another simple case. The 
parameter B(z) in Eq. ( |123| ) should be evaluated at the spiral tip position. As the spiral 
tip turns around the spiral core, B varies harmonically in time at the spiral rotation period 
and the spiral drifts. Since the direction of maximum excitability viewed from the spiral 
center is along the gradient of B, one concludes that the spiral drifts perpendicularly to the 
gradient of B, along an iso-excitability line. 



B. Drift in an external field and filament tension 



It has been reported in previous experimental and theoretical studies P4]J25[] 

that a spiral drifts when it is submitted to a constant external field. Interestingly, the spiral 
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was found to drift at a non zero angle with the applied external field. In presence of an 
external field E which couples to the activator u, the activator reaction-diffusion (|3|) becomes 

d t u = eV 2 u + f(u,v)/e-E.X7u (127) 

A simple way to determine the effect of E is to view the wave dynamics in a frame M 
which moves at velocity E. In such a frame, the supplementary gradient term in Eq. (|127|) 



disappears and u simply obeys the field- less Eq. (|]). However, the controller equation is 
modified. It reads, in the excited region, 

d t v = l/r e + E.Vv (128) 

The gradient term in Eq. (|128|) modifies the relation between the tangential tip velocity and 



the medium parameters. As shown below, one obtains instead of ( 37]) 





c t = c + c — — - + 7|| E\\ + 7 _l E ± (129) 

where En and E± are the external field component respectively parallel and orthogonal to 
the tangential tip velocity (measured in the frame M). Our sign convention is that E± > 
when it points toward the excited region of the spiral tip. The numerical coefficients 711 and 
7_i_ are determined below from a solvability condition, as we have now done several times. 
Before detailing this computation, we show that the spiral drift is a simple consequence of 
Eq. ( [1291) . As above, the wave tip motion is determined by Eqs. (|63|j6^) where now z = x + iy 



denotes the position of the wave tip in the frame M and q is given by (|129|) and depends on 



the angle between the instantaneous velocity (in the frame M) and the external field E. The 
form of the function Ri[ct] is a consequence of the front interface dynamics determined by 
(0) which applies in the frame M. Therefore, it still has the large core asymptotic form fl66|). 
Writing c t = c° + 5ce in ( |129| ) as a constant part c° independent of the external field and 
a small external field dependent part 5ce = J\\E\\ + J±E±, we can again copy the analysis 
of subsection ( |V A| ) and simply replace 5c q by 5ce- For a perturbed wave tip circle motion 
z = (R + eg(t)/c )e^ li+,/ ' (t) , this gives instead of © 
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2 2; C oR'i[ c t] 2 3 Cq-Rq $ C E 



q + ufq = ufdcE^^ = cuf- — ^ (130) 

e 2 e c - cl 

For definiteness, we suppose that the field E is parallel to the x-axis which gives, to lowest 
order in the perturbation, En = —Esm(uit) and E± = Ecos(uit). So, Eq. (|130|) is again 
found to be the equation of an harmonic oscillator forced at its natural frequency and the 
amplitude q of the oscillation diverges in time 

e 3 E 

q(t) = u)\t Yl\\ cos(u;it) + 7_l sin(a>ii)] (131) 



co-Ro 4 c - ct 

Comparing Eq. ( |131|) with the expression of q for a translated spiral (|125|) , one concludes 
that ( |131|) describes a spiral drifting away from the origin at constant velocity with 

3 E 

X ° = An 7 7|| 

4 Co — Ct 

3 E 

Vo = 7-1 Ro vit (132) 

4 c - ct 

The spiral drift angle Qr> with the external field is therefore 

tan(0 D ) = 7j./7[| (133) 

Several remarks can be made : 

i) formally, Or, is the angle between the drift velocity and the external field in the M frame. 
However, the drift velocity in the large core limit is dominantly produced by the time depen- 
dent variation of the spiral radius and is much larger than the velocity difference between the 
lab. frame and the M frame. Terms of the same order as the velocity difference between the 
two frames have been neglected in obtaining ( |131|) . It therefore makes no sense to correct 
9r> for this velocity difference. 

ii) A constant field produces a spiral drift because the r.h.s. of the components of the exter- 
nal field in the tip frame, E\\ = —Esm(uit) and Ej_ = Ecos(u>it), oscillate at the resonant 
frequency u)\. A sinusoidal external field oscillating at u e has components in the tip frame at 
u e + loi and io e — LO\. A spiral drift is therefore induced by an external field when it oscillates 



at twice the spiral frequency (u e = 2ui), as noted in previous studies [^2| 
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iii) As said previously, the derivation of ( |132|) breaks down when the spiral center is no 
longer close to the origin and the linearization giving ( |131| ) becomes illegitimate. This can 
be cured as said above, by introducing from the start the spiral center the motion of which is 
determined through the requirement that no secular terms appear on the r.h.s of Eq. ( |130| ). 

It remains to obtain ( |129|) and compute the parameters 711 and 7j_. We consider the spiral 
(in the M frame) in a Cartesian coordinate system attached to the wave tip as in Fig fj. 
As before, the front interface y/(x) simply obeys Eq. fll8| ) in the tip region. However, the 
controller concentration on the back interface is changed by the external field (Eq. ( |128| )) 
and this modifies the back equation ([T9|). 

We begin by computing the controller concentration on the back interface. The time 
dependence of the field components can be neglected since it is on the scale of the rotation 
period, Rq/cq which is much longer in the large core limit than the time scale of interest, 
the spiral width traversal time e/cg. Eq. ( |128| ) thus shows that in the excited region v obeys 

v(t, x - E ± t, y - E\\t) = v(0, x, y) + t/r e (134) 

The concentration Vb(x) on the back interface at the point (x,yb(x)) is related to the con- 
troller concentration vq on the front interface at the point (xf, Vf(xf) +c t t(x)) at a previous 
time t(x) with 

Xf = x — E±t(x) 
y f {x f ) + c t t(x) =y b {x) - E\\ (135) 

Xf and t(x) are functions of x, the considered point of the back interface which can be 
determined perturbatively for small external field. Writing Xf = x + 5x,t(x) = t (x) + 
St(x) , one obtains to(x) = [yb(x) — yf(x)]/c t ,Sx = —E±to(x) and 5t(x) = to(x)(— E\\ + 
E±dyf/dx\ x )/c t . Therefore the controller concentration at abscissa x on the back interface 
is equal to 

v b (x) = v - t(x)/r e = v + Vf{x) ~ Vb{x) [1 + (-£||/Q + E ± /c t ^\ x )} (136) 

c t r e dx 
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The last term is the modification of v on the back interface coming from the external field. 

When (|136|) is taken into account, the back equation in the tip region reads (using as 
before space variables scaled by e/co), 



d 2 y b 
dx 2 



[• • -}oid - B— [y f (x) - y b (x)\ 



-E^/ct + Ejct ^l\ x 



1 + &f 

ax 



1 3/2 



(137) 



where [• • denote the terms on the r.h.s. of Eq. (|2"2"|). When the front and back equations 
are linearized around the critical finger as yj(x) = Yf(x) + Syf(x),y b (x) = Y b (x) + Sy b (x) one 
obtains as before Sy/(x) = rjiSct/co (Eq. (|25|) and (p2|) ) and a modified equation for Syb(x), 



£b(5y b 



\old 



B c [Y f (x)-Y b (x)) 



-E\\/cq + E ± /c 



dx x 



1 + 



,dY h 



dx 



.3/2 



(138) 



Integrating both sides of Eq. ( 138 ), one obtains the solvability condition which replaces 
Eq. (M) 



Set 

Co 



[h + B c (-I 2 + I 3 )} -5Bh = B e (-EJco h + E ± /c 1 ± ) 



(139) 



where the constants I±, J 2 , /3 have previously been defined (Eq. (|36|)) and the new constant 
I± is given by the following integral 

3/2 

(140) 



+°° dYt 
dxi{x) [Y f (x)-Y b (x)]-f- 
o ax 



1 



'dY h 



dx 



8.431 



Eq. ( |139|) shows that (|129| ) holds with the following expressions for 7y and 7j_, 



7|| 



1± 



K 
BJ± 

Kh 



-1.177 



1.287 



;i4i) 



Changes of spiral core radius are the dominant effect in the large core limit and lead to a 
drift opposite to the field (711 < 0) as qualitatively argued in ( [p5|| ). We quantitatively find 
here that a counterclockwise rotating spiral drifts at an angle of about 132.5° with the field 
direction in good agreement with previous simulations [^5j as well as our own, as shown 
in Fig. [TJ| (the sign of 7j_ and of the drift angle would be opposite for a clockwise rotating 
spiral). 
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Finally, we note that the curvature induced motion of a weakly curved 3D scroll wave 
P6| , |3"T|| filament is directly related to spiral drift in an electric field. For a 3D filament 
(xq(s) , yo(s) , Zo(s)), we can choose a coordinate system with its third axis aligned with the 
filament tangent at s. Locally, the activator field can be written u(x — Xo(s),y — yo(s);t) 
with u(x,y;t) a two dimensional spiral wave. The two dimensional laplacian in Eq. ([|) 
acting on such a solution gives V| D tt — (x"d x + y"d y )u = V\ D u — kN.Vu where k is the 
filament curvature and iV the filament normal with kN directed toward the filament center 
of curvature. Therefore, enN acts as an external field E in the normal (x, y) plane. Since 
711 < and a spiral drifts opposite to the field direction, one concludes that curvature is 
destabilizing in the large core limit (negative line tension) and that a scroll ring grows. 
Moreover, it propagates normally to the plane of the ring at a velocity proportional to its 
expansion velocity since 7j_ > 0. The other laws governing filament motion can similarly 
be deduced by reducing the 3D dynamics to an effective 2D process. We defer, however, a 
detailed study of 3D dynamics in the large core limit to a future publication. 



VIII. MULTIARMED SPIRALS 

In this section we extend our analysis to the situation where several thin excited regions 
or 'spiral arms" rotate around a common core. Our main finding is that such multiarmed 
spiral waves are always linearly unstable in the large core limit. We confirm this finding by 
numerical simulation of the FitzHugh-Nagumo model for two-arm and three-arm spirals. A 



different conclusion has been reached in ref. [26] where multiarmed spiral waves were found 



by numerical simulation of the FN model, with a well-prepared initial condition, to be stable 
over windows of parameters in the large core limit. We shall comment at the end of this 
section, on the possible origin of this disagreement. 

Let us denote by qj{0) the coordinate of the tip of the j th spiral arm. We make the 
arbitrary choice that rotation is counter-clockwise and take the index j 6 [0, N — 1] to 
increase clockwise. The equation for the phases, ipj — @j ~ &it, are given by 
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di/ij/dt = -(e/coRo)ujiqj (j = 0, N - 1) 



(142) 



For simplicity, we consider an initial condition where the angular positions of the N spiral 
arms are uniformly distributed. To lowest order in e/(co-Ro), one can assume that the spiral 
arms rotate at constant angular velocity and that the phase difference between two successive 
arms remains constant: ipj = 2n/N. The equation that governs the motion of a given 

arm, say arm j, is essentially the same as the one governing the motion of a one-arm spiral, 
except that this arm interacts with the exponential recovery tail of the controller field v of 
arm j — 1, instead of its own recovery tail. Consequently, the equation of motion for arm j 
is simply obtained by replacing the interaction term mF(q(9) — q{9 — 271")) on the r.h.s. of 
Eq. |5| by rriNF(qj(9) — qj-i(0 — 2ir/N)), with m N defined in terms of the reduced period 
2tcRq/N. For a spiral with N arms, the wave tips are governed by the N coupled equations 

d 2 qj 



d6 2 



+ qj = m N F ( qj (9) - qj-x{e - 2tt/JV)) 



(j = 0,...,N-l) 



(143) 



where 



3/2 



m N 



3B c (bK) 
2(B C - Bff 



exp 



2ne 



bK 



3/2' 



c&Nt r \B c -B / 



(144) 



and F is the same function as for a one-arm spiral. 



A. Linear stability 



Let us first analyze the linear stability of an iV-arm spiral. Linearizing Eqs. |143| , we 
obtain 

^l + qj = a ( qj (9) - q^O - 2n/N)) (j = 0, N - 1) (145) 

where we have defined a = mjv-F"(0). The symmetry of the above system of linear equations 
implies that its solutions must be of the discrete Floquet-Bloch form 

qj = qexp(ik n j + Q n 9) (146) 
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where k n is the discrete Bloch wavector that takes on the values 



27TTI 

kn = ^ (n = 0,..,JV-l) (147) 



Substituting the above form into Eq. |145| , we obtain the eigenvalue equation 



nl + 1 = a 



2tt 

1 - exp ( - — (Q n + in) 



(n = 0,...,N- 1) (148) 



that determines the allowed values of Q n for each mode n and hence its stability. The two 
global translational modes, which are exact solutions of Eq. |148| for arbitrary a, correspond 



to Qi = —i and Qn-i = i. We restrict ourselves to considering the 2N — 2 other modes 
which correspond to the coupled translations of the individual spiral arms. The eigenvalues 
corresponding to these modes can be calculated perturbatively by expanding Q n in a power 
series in a about ±z. For brevity of notation, let us denote by Q„ the N—l eigenvalues 
obtained by expanding about Oi = +i for n = 0, 2, N — l, and by Q~ the ones obtained by 
expanding about fijv-i — —i for n — 0,1, N — 2. Substituting the power series expansions 

n n = ±z + a n± {1) + « 2 ^ (2) + ... (149) 

into Eq. |148| we obtain after simple algebraic steps 

Since the leading term in the expansion (|149|) is purely imaginary, the stability is determined 
by the sign of the real part of fi n (i). Eq. |150| implies that Re(0~Q\) > for n = or 
n > N/2 + 1, and Re(0^/ 1 x) > for n < N/2 — 1, and therefore that N-arm spirals are 
always unstable for N > 2. For the special case N = 2 and n = 0, Eq. |150| implies that 
Re(fi^ 1 - ) ) = 0, in which case the stability is determined by the sign of the real part of 
the next order term in the expansion, Re(f2^). The calculation at order a 2 yields that 
^q(2) = 7r/2=pi/2 and therefore that Re(f2^) = n/2 > 0. Thus the symmetric (n = 0) mode 
is always linearly unstable for a 2- arm spiral. In contrast, for the antisymmetric (n—l) 
mode, Qi = ±z remains solution for arbitrary a. We conclude that iV-arm spiral waves are 
always linearly unstable for iV > 1 in the large core limit. 
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The nature of the linearly unstable tip trajectories are simple to deduce from the above 
results. To be concrete, let us consider 2-arm and 3-arm spirals that we shall study in 
simulations below. For N = 2, aside from the two translational modes, there are two 
unstable modes corresponding to the complex conjugate pair 

ftjfc = vra 2 /2 ±i(l-a- a 2 /2) (N = 2) (151) 

Since this pair corresponds to n = 0, the two tips will move symmetrically (with equal radial 
displacements) about a fixed center of rotation. Furthermore, since the imaginary part of 
is slightly less than unity, the two tips will oscillate in and out of the unperturbed steady- 
state circle of rotation with a period slightly larger than the basic period T , and with an 
amplitude of oscillation that grows exponentially in time. For N = 3, there are four modes 
aside from the two global translational modes: a complex conjugate pair with a negative 
real part, which is stable, and the unstable complex conjugate pair 

fi± = VZa/A ± % (1 - 3a/4) (N = 3) (152) 



obtained by evaluating Eq. |150| for N = 3, where the tips move with equal radial displace- 



ments. As for N = 2, the finite imaginary part slightly smaller than unity implies that the 
tips will exhibit exponentially growing oscillations with a period slightly larger than To. 

In addition, a is typically much smaller than unity in the large core limit since the spiral 
period is large compared to the recovery time, tr <C Tq/N, and the spiral arms are only 
weakly coupled via the controller field v. Therefore, the instability of a multiarmed spiral 
should generically develop on a time scale much longer than T , especially for N = 2 since 
the real part of Qq scales as a 2 , instead as a for iV > 2. 

B. Numerical simulations 

In order to test the above predictions, we investigate numerically the stability of spiral 
waves with two and three arms in the FN model. We restrict ourselves to a range of 
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parameters where a one-arm spiral is linearly stable and rotates rigidly. We construct an 
initial condition for an iV-arm spiral, denoted by (u^, vn), by simply rotating N — 1 times 
by 2tt/N a one-arm spiral wave, which yields the expression 

JV-l 

u N (r, 0) = J> ( r > 9 ~ 2 ^/ N ) ~ ( N ~ l)«o (153) 

j=0 
N-l 

Mr, 0)= 5>(r,0- 2nj/N) - {N - l)v (154) 

j=0 

where (uq, Vq) are as before the resting values of u and v. Since the simulations are performed 
in Cartesian coordinates, and the edges have a negligible effect, each rotation of 27r /N 
is simply carried out by running the simulation of a one arm spiral for a time equal to 
Tq/N. The initial condition defined by (|153|) and ( |154j) deviates from the true steady-state 
solution of an iV-arm spiral by an amount proportional to v — v o on the wave fronts, which is 
exponentially small in the large core limit. Therefore, this initial condition can be considered 
as a slightly perturbed iV-arm spiral solution and is ideal for the present purposes. 

Results of the simulations are shown in Figs. [16] where we plot the normalized radial 
displacement of the wave tips, (rj(t) — R Q )/Ro, which corresponds to eqj/c in our analysis. 
We calculated the position of the iV wave tips by looking for the points of zero normal 
velocity along the spiral boundary defined by u = 0. This is equivalent to looking for the iV 
intersections of the curves u = and d t u = 0. We measured rj(t) from the instantaneous 
center defined by x(t) = Ef =1 Xj{t)/N and y(t) = I$ =1 yj(t)/N). All the main qualitative 
features predicted by our analysis are observed in the simulations. (We have not attempted a 
detailed quantitative comparison because our predictions are strictly valid outside the range 
of our simulations.) Firstly, during the initial instability, the center of rotation (x(t),y(t)) 
remains fixed in time and the radial displacements are equal for all tips. This implies that 
the symmetric n = is the most unstable one. Secondly, the radial displacements exhibit 
exponentially amplified oscillations, with the amplification rate depending sensitively on 
the steady-state period To and the number of arms, which both determine the parameter, 
a = F'(0)mN, entering in the predicted amplification rates (i.e. the real parts of fi^ in 
Eqs. |151| and [L52|) . In particular, Fig. [16] shows that the amplification is much slower in 
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(b) than (a), which agrees with the fact that To is about 1.46 times larger in (b) than (a). 
In addition, for the same parameters, the three-arm spiral in (c) is destabilized much faster 
than the two arm spiral in (b), in agreement with the fact that uin defined by Eq. |144] is 



larger for N = 3 than N = 2. Lastly, the period of the radial oscillation is slightly larger 
than T as predicted by our analysis. This can be seen for example in Fig. |iT)|(b) where the 
radial displacement of the tips exhibits 48 peaks over a time lapse of 50To. 

One interesting question is whether the instability of the symmetric mode saturates in a 
nonlinear regime. To explore this question, we have integrated Eq. |143| numerically for the 



symmetric mode by letting q\(t) = q 2 (t) = ...qiy(t) = q(t), in which case Eq. |143| reduces to 
a single equation for q(t). We investigated different values of N and m for the function F 



defined by Eq. |121| with a = 0.2. The results are shown in Fig. [17] for iV = 2 and N = 3, the 
plots for higher N being qualitatively identical to the plot for N = 3. These plots show that 
the bifurcation is subcritical. For all N > 2, the amplitude of oscillation increases linearly 
in time in the nonlinear regime. This comes about because in the forced harmonic oscillator 
equation for q the amplitude of the resonant forcing term F saturates when q becomes 
of order one. For small m^, averaging the forcing term over one period of the harmonic 
motion, gives the mean energy increase of the oscillator and accounts for the phenomenon. 
Interestingly, the cross-over from the linear to the nonlinear regime is qualitatively different 
for iV = 2 and N > 2. For N > 2, the slope of the envelope of the oscillations increases 
monotonously in time until it reaches a constant value in the nonlinear regime. Whereas for 
N = 2, the slope of the envelope increases non-monotonously with time. The FN simulation 
for N = 2 shows qualitatively the same non-monotonous increase of the envelope of radial 
oscillations with time as obtained by integrating the wave tip equation, as can be seen by 
comparing Fig. |16|(b) and Fig. |T7|(a)). This shows that even relatively fine details of the 
nonlinear instability of multiarmed spiral waves are captured by our analysis. In Fig. [16|(a), 
the oscillations grow too rapidly to their saturated values to observe this cross-over. 

One important consequence of the absence of a weakly nonlinear saturation of the un- 
stable symmetric mode is that the distance of closest approach between the wavetips (which 
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occurs at the minimum of each oscillation) decreases with time. The resulting highly non- 
linear regime is obviously not described by the wavetip equation ( |143| ) , which is only valid 
for small radial displacements of the wavetip compared to Rq. Results of the FN simulations 



show the complexity of the dynamics in this regime, as illustrated by Figs. [Lq and 19 



To conclude, let us contrast our results to those of Ref. [2~E\ where the stability of multiarm 
spiral waves was studied in a slightly different version of FN kinetics, but in a similar 
regime of weak excitability. When starting from sufficiently well-prepared initial conditions, 
multiarmed spiral waves were found to be stable when the period To was large enough to 
accommodate a finite number of arms around a single core. Moreover, it was observed that 
a spiral with N arms became unstable and decayed into a spiral with N — 1 arms when 
a transition line was crossed by decreasing T in the plane of T and the refractory period 
(defined as the minimum interval between waves in response to the lowest stimulus exciting 
the medium), with a separate line for each N. The main difference in our predictions is that 
steadily rotating multiarmed spiral waves are always linearly unstable for N > 2 for any 
parameters in this plane. Note however, that steadily rotating multiarmed spirals were not 
observed in |3B| when starting from randomly broken arms. 

We have actually checked that the instability predicted by our analysis, and observed in 
our FN simulations, also occurs in the FN kinetics studied in ||26|| . This is illustrated in Fig. 



20 for a two-arm spiral and k„ = 5.2, other parameters being chosen the same as in Ref. [26 



The main difficulty in observing this instability is that it develops extremely slowly when 
the spiral period is much larger than the refractory period, in which case defined by Eq. 
144| becomes exponentially small, and the time to observe the instability exponentially large, 



as a function of the ratio of the two periods. For example, for the parameter of Fig. |20|, the 
destabilization of the two-arm spiral already occurs over a timescale of about 10 rotations. 
For the value k g = 5 reported in Fig. 3 of [5f|, T is about twice larger than for k g = 5.2. 
Hence, the instability cannot be seen on a time scale of a few rotations. 
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IX. TOWARD SMALLER CORE RADII: A DISCUSSION 



We have seen in section [VI| that the large core equation of motion fl85|) leads to a mean- 
der onset frequency u>2 which is equal to half the basic spiral frequency, quite independently 
of the detailed form of the function F. It is interesting to identify the main subdominant 
effects which leads U2/0J1 to depart from 1/2 for smaller core radius (as shown in Fig. ID). 
The following two assertions underly the large core result: 

i) the tangential velocity and spiral tip rotation rate only depend on the instantaneous char- 
acteristics of the medium in which the spiral tip propagate (i.e. the relaxation of the tip 
velocity and rotation rate can be taken to be instantaneous), 

ii) the angular tip position is slaved to time (6 = uoit) i.e. the time interval between two 
successive passages of the spiral tip by the same angular position 9 can be taken to be 2it/ui 
and one can neglect the dependence of this time interval on the spiral path. 

A systematic discussion of corrections to the large core limit is beyond the scope of this arti- 
cle. We content ourselves here, in showing that corrections to i) and ii) both affect the value 
of U2/0J1 at onset. As discussed below, taking into account the non-instantaneous relaxation 
(i.e. corrections to i)) formally appear to give the dominant correction to the large core 
limit results. However, corrections to ii), although subdominant, seem the most important 
for the parameter range of Fig. |TD| and account semi- quantitatively for the numerical results. 

We begin by discussing i). The motion of the spiral tip is determined from the two 
relations (p3D,(p4|). The tangential tip velocity is determined by the dynamics of the close 
tip region ~ e/Ro, which is fast and independent of the spiral core size. The determination 
of the radius of curvature of the spiral tip trajectory involves however the dynamics of a 
whole intermediate region ~ (RoRf^) 1 ^ 3 and, as discussed in Appendix ||, this happens on 
a time scale td with u)\td/ ~ (Ru p / Ro) 1 ^ 3 ■ So one expects that this instantaneous radius 
of curvature, which we denote here by Ri to distinguish it from the steady state value Ri, 
adapts on a time scale td to changes in medium conditions. Short of solving Eq. ( |B2|) , a 
crude model of this effect is obtained by replacing the instantaneous Eq. fl64|) by 
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t d ^- + Rt = Ri[ct{{v})] (155) 



This gives instead of Eq. (ffoj) , the couple of equations 

d 2 q 



2 +q = 5R iCo /e (156) 



d9 

wiU + SRt = 5Ri (157) 

au 

where 6R4 is given by Eq. ( |7T| , |8^ ), as previously. At the linear level, Eq. ( |157| ) simply becomes 

uit d + 5R l = a -\q{9) - q(6 - 2tt)] (158) 



Searching for the eigenmodes of ( |156| , |158[ ) under the form q = Aexp{a6), gives the modified 
eigenvalue equation 

{a 2 + 1)(1 + auJxTd) = a[l- exp(-27ra] (159) 

The meander threshold is determined by requiring that Eq. (|159|) has purely imaginary 
roots a = iQ besides the two translation modes a — ±i. Perturbation around the large core 
{td — 0) result give the modification to the meander frequency at onset 

n = i - ^ (160) 

2 2n 

So relaxation effects lower the frequency ratio o^/^i below 1/2 and cannot account for the 
numerical observations reported in Fig.[T0|. 

In contrast, we show that improving on ii) lead to corrections in agreement with the 
numerical data. We parameterize the spiral tip position as in section ( |V A|) as z = {Rq + 
eq/co) exp[i{uit + ip(t))]. The angular tip position is 

9 = uJxt + ij;{t) (161) 

Beyond leading order, ip{t) is not negligible in ( |161| ) and the spiral period of rotation T 
depends on the spiral tip path. Eq. fl6"7|) gives ip = —uiq/R to dominant order (near 
the meander onset the other term in ( |67|) is of higher order, 5c q /c ~ {e/c R ) 5 ^ 3 using 
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Eq. (|70| ) |7lD and it can be neglected). This implies that it actually takes a time T = T + AT 
longer (shorter) than the period T of the steady spiral to return to the same 6 for outward 
(inward) displacements 

AT = To ^ f %(<j>) (162) 

CqHq JQ-2-k Z7T 

This reduces (increases) the interaction with the previous tip by ~ AT/r^exp(— Tq/tr) and 
causes the spiral tip to return sooner inside (outside) the core. This leads co^/^i to move 
away from 1/2 toward unity. In order to explicitly show this, we compute the variation of 
the tip trajectory radius of rotation 5Ri due to the spiral displacement taking ( |162|) into 
account. Comparing Eq. (|83|) and ( |162| ) gives 



SR 3B C f R c \ 
Ro 



2/3 



2bK [—) ex P(- T o/^)) 



F(q(6) - q(9 - 2k)) 



AT 2K + JB C 



TR B c 



(163) 



A modified version of fl85l) is obtained by substituting (|163|) in (|69|) 



q + ufq = mul \F[q(9) - q{6 - 2k)] - ?L 2K + JB ° \ (1 64 ) 
where the constant m is given by ( jS7D and AT depends on the tip trajectory (Eq. ( |162| ). 



The linear version of (|164j) is (where we can replace time by angular position), 

^| + q = a (q(6) - q(6 - 2k)) - /3 f ^ q{cj>) (165) 

where a = mF'(0) as before and (3 = m(To/TR)(e/coRo)(2K + 2JB C )/B C . The eigenmodes 
of ( |165|) are of the form q(9) = Aexp(ad) where a is a solution of 



cr 2 + l = ( a _-^_)(i_exp(-27ro-)) (166) 
2na 

The meander onset corresponds to the critical value a c where ( |166| ) has purely imaginary 
roots a = iQ (besides the two translation modes a = ±z). For small (3, first-order perturba- 
tion around the j3 = values gives 

Qc = 8 " 3^ (167) 

2 6K Z 
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This shows that the correction term in ( |164f ) lowers the threshold for the meander instability 
(i.e. plays a destabilizing role) . More importantly, it increases the frequency ratio Q = u 2 /uJi 
at meander onset, as announced. 

The frequency shift predicted by ( |168| ) can be compared to the numerical results of 
Fig. [It]. Using the lowest order threshold estimate mF'(O) = 3/8, one obtains Q — 1/2 ~ 
t / (coCtTnF' (0)) . With the estimate F'(0) ~ — .581n(A), the frequency shift is found to be 
of the same order of magnitude as the one measured. This semi-quantitative agreement 



leads us to think that, for the parameters of Fig. |10|, the correction (|168|) is the main 



effect and that the correction (|160|) is still numerically smaller than (|168|) . Of course, for 
spirals of sufficiently large core this should cease to be true, the correction ( |160| ) should 
become dominant and co^A^i is expected to drop below 1/2 before ultimately reaching its 
asymptotic value. Unfortunately, a numerical check of this non-monotonic behavior would 
require simulating spirals of very large core radius. This appears a difficult task with present- 
day computers. 



X. CONCLUSION 

We have developed an analytical approach to spiral waves close to the line OR where 
the spiral rotates around a large core and in the free boundary limit where the medium 
exhibits an abrupt response to a stimulus (e <C 1). The main ingredient of our analysis 
has been to note that in this limit the entire wave tip can be treated as an essentially rigid 
body, the slow motion of which is controlled by the local spatial gradient of excitability 
in the medium in a way that can be precisely deduced from the starting reaction-diffusion 
equations. This has provided a simple understanding of the spiral tip motion and a precise 
reduction of its dynamics to that of a single point. This has allowed us to describe the 
Hopf bifurcation nature of the meander instability and to derive simply, but with precise 
asymptotic estimates, spiral drift due to spatial or temporal variation of excitability, or due 
an imposed external field. This last computation determines in particular the drift angle of 
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the spiral with the external field and also the parameters governing the motion of an average 
scroll wave filament (curvature has been found to be destabilizing in the large core regime). 
In addition, our analysis has allowed us to elucidate a generic instability of multiarmed spiral 
waves that was previously missed in numerical simulations in the large core limit because it 
develops very slowly. 

The present analysis can be compared with several previous analytical approaches which 
have provided insights into spiral wave dynamics. As already noted, a phenomenological 
kinematical model of spiral wave dynamics [|l8j has been proposed several years ago and 
has succeeded in capturing many aspects of spiral wave motion. It differs from the present 
approach not only because its parameters need to be adjusted and cannot be obtained from 
the underlying reaction-diffusion equation but also more fundamentally because, here the 
dynamics of the spiral tip is reduced to an ordinary differential equation (ODE) and drives 
the motion of the rest of the spiral arm whereas in the kinematical model of |18| the tip 



motion follows from that of the whole curve. Moreover, the tip motion is described here 
in a different way, by the tip rotation rate and not by a growing or retracting velocity 
18 1 . Another notable approach is based on normal forms [|16| . As in our case, the 



as m 



tip motion is described by ODE. The normal from approach postulates the existence of 
a Hopf bifurcation and it describes its coupling to the spiral translation modes and the 
resulting tip motion based on general symmetry arguments, close to the resonant case where 
the meander frequency u 2 is equal to the basic spiral rotation frequency lo\. The present 
approach is restricted to a particular limit but makes more specific predictions. Besides 
providing determined parameters in the reduced equation which gives, for instance, the drift 
angle with an external field, it has the advantage, in our view, to provide an understanding 
of the physical mechanisms responsible for the very existence of spiral waves and of their 
dynamics, be it meander or drift due to external action. 

Extensions of the present work can be considered in several directions: 
- It would be interesting to extend the analysis to slightly more excitable media to capture 
hypermeandering or at least the change from inward to outward petals (i.e. the line U\ = u^)- 
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This would require going beyond our adiabatic approximation and considering the dynamics 
of the intermediate region. 

- The large core nature of the spiral rotation (i.e. the proximity of the line dR) is an essential 
element of our approach but several of our arguments do not really require sharp front and 
back interfaces (i.e. e C 1). This is certainly true for the —3/2 divergence of the spiral 
radius divergence near the line dR which only requires that the spiral normal velocity and 
curvature be related by an eikonal equation of the type of Eq. (|) on a sufficiently large 
scale. This is also the case for the validity of the adiabatic approximation. Thus, it appears 
that a computation of critical fingers and of the allied solvability conditions at finite e would 
provide an extension of our reduced description to the neighborhood of the full line dR. 
This would not accurately describe meander(since dR and dM are close only for e< 1) but 
would allow a simple quantitative description of spiral drift and other phenomena along this 
line. 

- Finally, it appears possible to extend some of our calculations to scroll waves in 3D as 
succinctly described for filament motion in section ( |V11 Bp . Hopefully, this will not only 
provide definite coefficients in the average filament equations of motion, but it will also 
provide a better understanding of the dynamics and instabilities of 3d scroll filaments pT7| , f4"8 
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APPENDIX A: THE TIP BOUNDARY LAYER 

It has been noted in section ( [IV A| ) that the solutions of the free boundary problem (^^) 
are continuous as well as their first two derivatives but have a discontinuous third derivative 
at their tip. We show in this appendix that this weak non-analyticity can be taken care of 
by introducing a boundary layer of size e/^/co ~ e 5 ^ 6 near the wave tip (i.e. smaller than 
the tip radius of curvature of size e/co ~ e 2 / 3 and larger than the interface width ~ e). 

We restrict ourselves to analyzing the case of a critical finger. We take the interface 
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width e as length unit. We first consider the sharp interface case. We find it convenient to 
parameterize the interface as x = h(y) instead of y = yj/^x) as in the main part of this 
article. In the vicinity of the wave tip (h'(y) <C 1), Eq. (|5|) reduces to 

c h\y) = c(v) - h"{y) (Al) 

The non analyticity of the interface is a direct consequence of the non-analyticity of c(v), 

c(v) = c = -a(v -v s ), y>0 

OLE 

c(v) = c - 2y + • • • , y < (A2) 

CQTe 

gives 

Ky) = -[\{yco?-\{ycvf + •••], y>0 

Co 2 o 

Kv) = -ikvco) 2 - ~(1 - 2B c )(yc ) 3 + •■■], y < (A3) 
Co 2 o 

When one takes into account the finite width of the interface c(v) becomes a rapidly but 
smoothly varying function in the wave tip neighborhood. For a shape moving at Co along 
the ^/-direction the two reaction-diffusion equations (|3|||) become 

- c c\w = V 2 m + f(u, v) (A4) 
-c d y v = eg(u,v) (A5) 

For notational simplicity, we consider functions / and g of the form f(u,v) = F(u) — 
v, g(u,v) = u — 7] with the stall concentration v s = and the corresponding rest state at 
u = 0,v = 0. We choose F(u) = —Au(u — l)(u — 2) for illustrative purposes which gives 
A = 2Ai]. To study the tip neighborhood, it is convenient to use instead of x, the displaced 
coordinate z = x — h(y) where x = h(y) is a line in the interface transition region (for 
definiteness, one can take an iso-w line, for instance the line u = 1 with the above choice of 
F). Eq. QAg ) then reads 

- ca[d y - h\y)d z ]u = d 2 z u + [d y - h'(y)d z fu + F(u) - v (A6) 
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Eq.(Al,A2 



The controller field v is assumed to remain close to the stall concentration and in the tip 
neighborhood h! is small. Thus, at dominant order, Eq. (|A6|) reduces to 



d 2 z u + F(u) = (A7) 

which has a standing front solution u^(z) which goes from = at z = — oo to = 
at z = +00 ( u^(x) = 1 + tanh(xy^4/2) for the above choice of F). At next order, one 
obtains, 

d 2 z u (1) + F'(u (0) )u^ = (coti(y) + h"(y))d z u® - h'\y)d 2 z u^ + v (A8) 

Integrating both members of Eq. (|A8|) with the zero- mode d x u^ gives the solvability condi- 
tion, 

c h'(y) = -h"(y)-c(y) (A9) 

with 

Jdz vd z u^ 

C{V) = Jdz (A10) 

When v has negligible variations in the interface width, Eq. ( |A1U| ) gives back the sharp 
interface result with c(v ) = — av and 

Q = (Jdz {d z u^)f (AU) 

On the contrary, in the tip region, v varies in the interface transition region and the integral 
term in ( |A10|) needs to be more carefully evaluated. To lowest order, the field v on the 



interface is obtained by integrating Eq ,( [A5| ) 

v(z,y)=v + — dyi u (z + h(y) - h(y x )) (A12) 
Co Jy 

When ( |A12| ) is substituted in ( [A10| ), one obtains a smooth function c(y) 

c(y) = -av = / d yi T(h(y) - h{ Vl )) (A13) 

Cn Jy 
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with 

t( W ) = r dz d ^ u j z+w) (am) 

J-oo « M 

Eq. ( |A14| ) gives a smoothly varying function (for instance with the above choice of F, 
T(w) = 1/2 [expOyC4/2)/ smh(w^A/2) - w^A/2/ sinh 2 (w^4/2)]) instead of the Heavy- 
side function of the sharp interface limit. To make further progress, we assume (and check 
afterwards) that h(y) = coy 2 /2 + f](y) with 77 a small correction in a neighborhood of the 
spiral tip that can be neglected in evaluating the integral term in Eq. (|A13|) , 



r+00 f+00 I 

/ d yi T(h(y) - h{y x )) ~ / d yi T(c y 2 /2 - c y 2 /2) = —S{y^) (A15) 

Jy Jy ^/Cq 

Eq. ( |A9| ) then gives for the tip profile correction 

c 2 y + cov'(y) = - t -^n L s (vVco) - v"(y) ( A16 ) 

c 



Comparing the different terms, one obtains that a consistent scaling is y ~ 1/ ^/co and 
V ~ v^o wn i cn gi ye the size of the boundary layer (note that, here, our unit of length is 
the interface width e) and the magnitude of the shape correction in the boundary layer. 
This legitimates the neglect of 77 in ( |A15[ ). In the scaled variables, y = Y ' j ' ^/cq and rj(y) = 
fcQH{y«Jco), the equation for the tip profile correction is 

d 2 H 



dY 2 



B c S(Y) + Y = (A17) 



where the function S(Y) is defined by (|A15|) and ( |A14|) (the second term on the l.h.s of ( A16|) 



is of higher order). The behaviors of S at infinity, S(Y) — > at Y = +00 and S(Y) — > —2Y 
at Y = —00, give the corresponding asymptotic behaviors of H, H(Y) — > —Y 3 /6 at Y — +00 
and H(Y) — > — Y 3 (l — 2B c )/6 at Y — —00. These precisely match the different small y 
behaviors (Eq. (|A3|) ) of the sharp interface description. It shows that H{Y) interpolates 
smoothly between these different behaviors. 
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APPENDIX B: DYNAMICS OF THE INTERMEDIATE REGION 

The analysis of spiral dynamics that we have developed in the main part of this arti- 
cle makes a crucial use of an 'adiabatic' assumption. Namely, that for changes of medium 
parameter on the time scale of the spiral rotation period T ~ Rq/cq, the instantaneous 
motion of the spiral tip can be taken to be that of a spiral tip moving in an steady medium 
with characteristics invariant in time and identical to those of the changing medium at the 
considered time. In the large core limit, the spiral is described by matching three regions : 

- a close tip region on the scale of the tip radius R tip = e/c which determines the spiral tip 
tangential velocity, 

- an intermediate region of size {RqR^) 1 ^ which determines the instantaneous radius of 
curvature of the tip trajectory 

- and finally, an outer scale the dynamics of which is driven by the previous two regions. 
The close tip region relaxes on a time scale which is independent of the spiral radius and 
which therefore clearly becomes short compared to the spiral period T for a spiral of suf- 
ficiently large radius. In this appendix, we show that the intermediate region relaxes on a 
time scale T (Rup/Ro) 1 ^ 3 which is also much shorter than the rotation period T in the large 
radius limit. This justifies our adiabatic assumption. 

We first write the dynamic equivalent of the static BCF equation ([H]), that is the motion 
of a curve governed by (||) using polar coordinates 



As for the static case, it is convenient in the intermediate region to introduce the rescaled 
variables y and £ with Of = uj]t + ye/(c R Q ) and r = R + [2_R (e/c ) 2 ] 1//3 £. Expanding the 




(Bl) 



square root in ( |B"T| ) and keeping terms of the dominant order gives 




(B2) 



with 
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2 -l/3 



/ Cq-Rc 



2/3 



Q — c 



(B3) 



V e J c 

Eq. ( P2|) is the dynamic equivalent of the static equation fl43|) determining the shape of the 
intermediate region. It shows that the characteristic time to adapt to changes of a (e.g. of 
c t ) for £ of order unity (e.g. for the intermediate region) is Rq/cq [e/ {cqRq)} 1 ^ . It is shorter 
by the factor (R tip / Rq) 1 ^ 3 than the rotation period as announced above. It is however worth 
pointing out than this is larger than the time one may have guessed, namely the length of the 
intermediate region divided by the velocity Cq. The reason is that in the intermediate region 
the interface is almost radial. As a consequence, the advection velocity is much smaller than 
Co and advective effects become comparable to diffusion-like effects due to surface tension 
(i.e. the last two terms in (|B2|) are of the same magnitude). 



APPENDIX C: SOLVABILITY INTEGRALS AND FUNCTIONS : SOME 

RELATIONS 



In this appendix, we recapitulate the definitions and give additional information on the 
several functions and integrals which have been introduced in the evaluation of solvability 
conditions. 

The linear operators considered are £ / and Cb which comes from the linearization of the 
front and back equations around the critical finger in the tip region 



with, 



C 



f 



dx 2 



-2 + 3 



1 + 



dY, 



f 



dx 



C b = _ a { x )~J~ ~ K x ) 
dx z dx 



a(x) 



2 + 3[l-£? c (Y f (x)-Y b (x))} 



b(x) =B C [1 + 



JY h 



6\2 



3/2 



dx 



1/2 



dY f d 
dx dx 



dx 



1/2' 



dY 
dx 



(CI) 

(C2) 
(C3) 



(C4) 
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Yf(x) and Y{,(x) are the critical finger front and back interfaces which satisfy (|T8|) and 
(|T9|) with B = B c = 0.5353 ■ • -. The small x behaviors of these different functions are 
Y f (x) = V2x+x/3+- ■ Y b (x) = -V2x+x(l-2B c )/3+- ■ •, a(x) = -3/ \2x) + y/2B c / \/x+- ■ ■ 
and b(x) = B c /{2xf/ 2 + ■■■. 

The zero-mode £(x), x > 0, of the adjoint of L b is the solution of 



4(0 = f| + ^-(a(x)O - b(x)t = 



(C5) 



which tends to when x — > +oo. It is here normalized by imposing the supplementary 
condition sup 2 . >0 [^(x)] = 1. A local analysis determines the behavior of £(x) for small x, 
£(x) = £'(0) [x — B c /^/2x 3 ^ 2 ln(x) + •••]. Eq.(p5|) has been solved numerically by a finite- 
difference scheme on a non-uniform grid (with a step size decreasing to zero at small x). A 
graph of the obtained solution is shown in Fig. ^. The computed value of the derivative of 
£ at the origin is £'(0) — 4.441. An exact relation between £'(0) and a weighted integral of 
£ is obtained by integrating ( |C5|) between x = and x = +oo, 



£'(0)/2 



dx £(x)b(x) 



BJ 



c J 4 



(C6) 



The verification of ( |C6| ) serves as a check of our numerical computation. 

To evaluate the solvability conditions, besides Yf, Yb and £, the solutions of the following 
inhomogeneous equations with the linear operator £/ (Eq. ( |Cl| ))are needed, 

1 , W\ a 



i + 



dx 



tji(0) = 0, ^(0) = 1/3 

3/2 



= (y>(x) - n(x)) 



772(0) = 0, rj 2 (x) ~ for x 1 

3/2 



dx 



^,o(0) = 0,^(0) = 2/3 



dx 



'dY, 



f 



dx 



77e(0) = 0, r/ e (x) ~ y / x/2 for x < 1 



(C7) 
(C8) 

(C9) 
(C10) 



They are plotted in Fig. |lq . 

For the evaluation of the different solvability conditions, it is useful to compute the 
following integrals, 
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h 



+00 



+00 



+00 



dx£(x 

dx£(x 

dx£(x 

dx£(x 

dx£(x 
dx£(x 



L v,0 



+00 



dx£(x 



+00 



dx£(x 



dx£(x 



\ dx J 



Tjl(x 



1 



'dY h 



2.771 

3/2 



dx 

[Y f (x) - Y b {x)\ 



3.814 



1 + 



dx 



3/2 



~ 7.708 



14 

7/2(2 
dY b 



dx 



3/2 



dx 



4.1476 

3/2 



dx 



Vv,0 



1 + 



v dx 
dx 



3/2 



6.306 



-2.118 



12.553 



[r,(x) - n(x)] 



di2 

dx 



1 + 



dn 
dx 



3/2 



dn 

dx 



4.476 



3/2 



8.431 



(Gil) 



Exact relations between some of these integrals can be obtained by using symmetry 
transformations of known action on the interfaces. For instance, under dilation the critical 
finger front and back become Yf a = Yf(ax)/a, Y b ^ a = Yf ) (ax)/a and obey scaled versions of 
Eq.QTSp and 



d 2 Y f , a 
dx 2 



d 2 Y h 



b,a 



dx 2 



a 



a 



1 + 



1 + 



,dYi 



f,a\2 



dx 



1 + 



,dY, 



f,a\2 



dx 



3/2' 



AY h 



b,a \2 



dx 



[l-B c a[Y f , a (x)-Y b , a (x)}} 



1 + 



JY. 



b,a \2 



dx 



3/2 ' 



(C12) 



(C13) 



An expansion around a = 1 gives Yf >a = Yf + (a — l)5Yf + - ■ ■ where 5Yf(x) = xYl(x)—Yf(x). 
Similarly, one has Y b;CX = Y b + (al)SY b + ■ ■ ■ with 5Y b (x) = xY b '(x) — Y b (x). Expanding 
( P12| , |C13|) in the same limit shows that SYf and SY b obey the following linear equations , 



C f (SY f ) 



AY, 



/^2 



dx 



1 + 



AY, 



/^2 



3/2 



dx 



(C14) 
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C b (6Y b ) 



1 + 



,dY 



&\2 



dx 



[1-2B C [Y f (x)-Y b (x)]-B c SY f ] 



,dY 



&\2 



dx 



3/2 



(C15) 



Eq.( |C14|) shows that SYf = r/i — r] 2 ( |C7| , p8|) since SYf = (as can be checked from its 
explicit expression). Then, multiplying both sides of (|C15|) by and integrating from 
x = to +oo gives the searched relations between the above integrals 



h + h- B c (I 2 + 2/ 3 - h) = 



(C16) 



Using rotational symmetry in a similar manner, one obtains that rj e (x) = x+YfdYf / dx—1 
and the other relation 



h + B C (I ± - I e ) = 



(C17) 



In the analysis of meandering, there appear several functions of the tip displacement: 

3/2 



1 



dY t 



f 



dx 



Vv,d{°) = ( C18 ) 



I 3d = dx£(x) [Y f (x + d) - Y b (x + d)} Q(x + d) 
Jo 



1 + 



'dY, 
dx 



3/2 



r+00 / 
Iv,d = J o dx^(x) T] Vjd 1 + [ 



'dYb 
dx 



3/2 



(C19) 



For large \d\, these functions tend toward constant values, toward when d — > —00 and 
I34 2/4/ 'B c and I v>( j — > 2I 5 /B C when d — > +00 Their behavior for small displacements 
(\d\ 1) is nonanalytic. From the small x behaviors of Yf,Y b and £, one obtains, for 
< \d\ < 1, 

r(o) 
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-d\n(\d\) + 0(d) 



r}vA x ) = Vv,o(x) - d\n(\d\) + 0(d) 



(C20) 
(C21) 



The expansion (|C21| ) of rj V) d for small d gives for I, 



v,d 



L d = L -hdln(\d\) + O(d) 



(C22) 
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The function F which measures the spiral self-interaction, F(d) = (I^,d + B C I V ^)/I^ — J, has 
therefore a singular expansion for 0< d< 1, 

F(d) ~ - ^ dln(\d\) ~ -.576 dln(|d|) (C23) 

(where we have used ( |U6|) which shows that the singular contributions of I^d and B C I V ^ are 
equal). For d — > — oo tends toward —J ~ —1.872 and for d — > +oo approaches 
2(h + I 4 /B c )/I 3 -J~1.774. 

The singular behavior of F(d) at small d disappears when the controller field diffuses. 
For small diffusion (Yf(x + d) — Y b (x + d))Q(x + d) is simply replaced in (|C18| |C19| ) by the 
smoother function YS(x + d) 

YS(x;£ D ) = / expC-^-^-JCF/CaO -Y h {x')) (C24) 

For ffl < 1, one can check that the singular behavior of the self-interaction function is 
cut-off at d ~ to and the small distance behavior of F(d; to) is 

F(d; t D ) ~ -.576 d \n(t D ) (C25) 

In the other limit to 3> 1, for a; of order one, YS(x;to) — 1/B C + 2x/(^B c t D )- The 
corresponding behavior of h,d-i D and I v ,d;i D at small d is h,d-e D = I%,Q;i D + 2I±d/ '(B c ^/nto) + 
• • •> Iv,d-,i D = Iv,o;£ D + 2I 5 d/ (Bcy/wto) + • • •■ So the small distance behavior of the spiral 
self-interaction function is , for I ^> 1 (but still smaller than the scale of the matching 
region) , 

F(d; to) ~ 2 h u + ^ ^- ~ 2.06d/^ (C26) 

B c ^7ll 3 t D 

The derivative at d = of the spiral self interaction function has been computed numerically 
for intermediate values of to- It is plotted in Fig. |9|. 

APPENDIX D: NUMERICAL INTEGRATION OF THE WAVE TIP EQUATION 



In this appendix, we describe a simple scheme to integrate numerically the equation 
of motion for the wave tip (j85|) , which is convenient to rewrite as a system of first order 
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ordinary differential equations 

de =p (D1) 

^ = _g + mi ?( g (0)_g(0_ 27 r)) (D2) 

The difficulty of integrating this equation comes from the fact that the translational invari- 
ance of the underlying reaction-diffusion equations remains present in Eq. flSSj), and hence 
in Eqs. (|Dl|) - (p2|) , which are invariant under the transformation 

q(6) = q(6) + Ae ie + c.c. (D3) 

where A is an arbitrary complex amplitude. It is therefore desirable to develop a numerical 
scheme that discretely preserves this symmetry in order to avoid spurious discrete effects 
resulting from the coupling of this translational mode to other modes. To see how to 
construct such a scheme, let us first consider the case the second term on the r.h.s. of 
Eq. ( P2|) is absent. In this case these equations describe simple harmonic motion with a 
constant energy ~ \A\ 2 . It is well-known (and simple to show) that the simple Euler explicit 
scheme 

q n +i = q n + hp n (D4) 
Pn+i =Pn~ hq n (D5) 

where h is the time (angle) step, does not conserve energy, but rather pumps energy into 
the motion. As a result, it leads to an unbounded increase in A at long time. Therefore, 
this scheme violates in an obvious way the symmetry that we would like to preserve here. 
In contrast, the modified (Euler-Cromer) scheme 

Pn+i =Vn~ hq n (D6) 
Qn+i = Qn + hp n+1 (D7) 

exactly conserves the energy of harmonic motion. This can be seen by substituting the 
ansatz p n = p r n and q n = q^r 11 into Eqs. ( |D6| ) and ( |D7| ). Nontrivial solutions then exist 
only if 
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r = 1 - hr/2 ±ihy/l- h 2 /4 = e ±l£Xti (D8) 

where 



, hJl-h 2 /4\ 

Ae ^ (-t^jt) < D9) 

Hence, q* = Ae inAe + c.c. is a solution of Eqs. ( |D6| - |i)7| ) with constant A. Let us now extend 
this scheme to the case where the second term on the r.h.s. of Eq. ( |D2| ) is included, by 
simply letting 

Pn+i =Pn~ hq n + hmF(q n - q n - N ) (D10) 
Qn+i = Qn + hp n+ i (Dll) 

Now the key point is that g* = Ae mAd + c.c. remains an exact solution of these equations 
only if q* — q*_ N = 0, and thus 

A6 = 2ir/N (D12) 



This condition together with Eq. ( P9|) then uniquely fixes the step h for a given number of 
time steps, N, per basic period of 2ir. After simple algebraic manipulations, we find that h 
should be equal to 

h = 2 sin(7r/iV) (D13) 



In summary, our integration scheme is uniquely defined by Eqs. ( pi(J| ) and ( [01 1| ) with 



h given by Eqs. (|D13|) . For an arbitrary value of N, this scheme is invariant under the 
transformation 



q n = q n + Ae mA0 + c.c, (D14) 

which is the direct discrete analog of Eq. (|E>3|) . A solution of a desired numerical accuracy 
can then be obtained by choosing iV sufficiently large. 
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FIG. 1. Plots of the propagation (dP), rotor (dR), and meander (dM) boundaries in the 
parameter space e (the ratio of the fast activator to the slower controller time scale) and A (the 
medium 'excitability' defined as v s — vo) for the numerically simulated FitzHugh-Nagumo kinetics 
(/(u, v ) = 3u — u 3 — v, g(u, v) =u — 8). Our analysis predicts that the three boundaries approach 

smoothly the origin without crossing as e — ► 0, with A p e^lne for dP, A c ~ e 1 / 3 for dR, and 

A m - A c e 5 / 9 /(hie) 2 / 3 for dM. The prediction A c = (2 1 / 2 4e/S c ) 1 / 3 with B c = 0.535 for the 

dR boundary is in good agreement with the simulations. 
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FIG. 2. Surface plots of u and wave tip trajectories (thick solid line) illustrating in (a) a 
retracting wave for 6 = —1.4 and e = 0.27, in between the dP and dR boundaries, and in (b) a 
large core meandering spiral wave for S = —1.4 and e = 0.18, close to the dM boundary. 
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FIG. 3. Sketch of the spiral tip region and unsteady tip trajectory (solid lines). (r,9) denotes 
the polar coordinates of the wave tip with respect to the fixed steady-state center of rotation O. 
Ri is the instantaneous radius of curvature of the unsteady tip trajectory about the instantaneous 
center of rotation Oi, with Oi = O and Ri = Rq for steady rotation, and Ru p is the radius of 
curvature of the boundary between the excited and recovery regions of the medium at the tip. c t 
denotes the instantaneous tangential velocity of the wave tip along this trajectory, with q = oj\Rq 
for steady-state rotation. The coordinates 5r = r — Ro = eq/co and tp = 9 — uj\t measure the 
radial and angular departure from steady-state rotation, respectively. The cartesian coordinate 
system (x,y) that moves with the wave tip is also shown with the y-axis parallel to c t . Uf(x) and 
Ub(x) denote the instantaneous wave-front and wave-back boundaries. Finally, d = q(9) — q(9 — 2ir) 
measures the radial displacement of the wave tip after one 2tt rotation. 
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FIG. 7. Schematic plot illustrating: (a) the variation, of the controller field v on the instanta- 
neous wavefront, Vf (solid line), resulting from the previous passage of the spiral wave at the same 
angular position with the tip displaced radially outwards by d. The dashed line in (a) indicates the 
variation of v on the waveback, Vb, at the time of the previous passage of the spiral. The solid and 
dashed line in (b) represent the spiral boundary at the present time (solid line) and at its previous 
passage (dashed line). Note that the excitability averaged along the instantaneous wavefront is 
higher than for steady-state rotation due to the radial displacement after one rotation. Our for- 
malism provides a rigorous procedure for calculating how the instantaneous tangential velocity of 
the wave tip changes in response to this spatially varying excitability. 
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FIG. 8. Graph of the of the spiral self-interaction function F(d) vs. tip displacement d for 
different diffusion lengths : in = (solid line), in = 1 (long-dashed line) and Id = 3 (short-dashed 
line) . 



87 



3.0 



2.0 

F(0) 



1.0 



0.0 




0.0 5.0 10.0 

FIG. 9. Derivative at d = of the spiral self-interaction function F(d;£r>) vs. the diffusion 
length ip. The dashed line shows the large Id (Eq. |C26| ) asymptotic behavior. 
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FIG. 10. Plot of u^/^i vs co/eRo obtained by simulations of FitzHugh-Nagumo kinetics with 
f(u,v) = 3u — u 3 — v and g(u,v) = u — 5 (solid line and circles). Co = (S 3 — 3(5) /2 1 / 2 . The dash 
line represents the extrapolation to the asymptotic limit u^/wi = 1/2 predicted by our analysis. 
These simulations were carried out using a second order accurate direction implicit scheme with 
dx/e = 0.33 and dt/e = 0.1. The insert shows an example of a large core meander pattern for 
e = 0.180 and 5 = -1.4, where u 2 /uji ~ 0.67. 
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FIG. 11. Plot of q{9) — q{6 — 2tt) vs 9/2ir obtained by numerical integration of the wave tip 
equation with F defined by Eq. |121| and a = 0.2; (a) m = 0.42, and (b) m = 0.6. The onset of 
meander for this function F corresponds to m c = 0.3902. 
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FIG. 12. Plot showing saturated oscillations of q(6) — q{9 — 2ir) vs 6/2ir for different m. 
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(a) 



(b) 



FIG. 13. Comparison of large core meander trajectories obtained: (a) by plotting a two-radius 
epi-cycle (solid line) with O = 0.735 and pi/po = 1/5 and the predicted three-radius epicycle 
(dashed line) with Q = 0.735, pi/po = 1/5, and P2/P1 = (1 — ^)/(l + and (b) by simulation of 
the FN model with e = 0.18 and 5 = —1.4. The value of Q and the ratio pi/po used as input in 
(a) were extracted from the simulation in (b). The total time in (a) and (b) is about 3Tq. 
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FIG. 14. Simulation of the FN model of Fig. |T| (e = 0.185, <5 = -1.41) with an external field 
added as in Eq. ([127D with E = 1.0 10~ 3 . The wave tip trajectory is shown (bold line) as well as 
surface plots of u showing the spiral position at the end of the simulation. The spiral is found to 
drift at about 135° with the field in good agreement with the asymptotic prediction of 132, 5° 




FIG. 15. Graph of the functions rj\(x) (solid line), rj2(x) (dotted line) and rj v fi(x) (dashed line) 



defined in Eq. ( |C7 
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FIG. 16. Plots of the radial displacement of the wave tips vs time for two-arm (a and b) and 
three-arm (c) spirals. The highly nonlinear symmetric meander dynamics of three-arm spirals (Fig. 
|j~8|(b)) is destabilized at large enough time leading to the elimination of one arm at boundaries. In 
contrast, the symmetric meander dynamics of two-arms spirals is stable on the time scale of our 



simulations despite the collisions illustrated in Fig. 19 
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FIG. 18. Simulations of the FN model showing the wave tip trajectories during the initial 
development of the instability of multi-arm spirals for: (a) a two-arm spiral with e = 0.445 and 
5 = —1.24, and (b) a three-arm spiral with e = 0.445 and 5 = —1.25. Each line type (solid, dashed, 
or long-dashed) corresponds to a different wave tip trajectory. 
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FIG. 19. Sequence of surface plots of u and superimposed wave tip trajectories (thick solid 
lines) illustrating the highly nonlinear collision of wave fronts occurring during a symmetric mean- 
der pattern of a two-arm spiral. Simulation parameters are 5 = —1.24 and e = 0.445. Frames (b), 
(c), and (d) are at t/e = 20, 25, and 40, respectively, with t measured from frame a. Note that an 
exchange of wave tips and spiral arms occurs during the collision in (c), such that the wave tip of 
the downward moving arm in (b) is at the end of the upward moving arm in (d) and vice versa. 
This exchange produces the sharp pivot turns around the small inward meander petals. 
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FIG. 20. Plot of the radial displacement of one of the spiral tip vs time showing the instability 



of a two- arm spiral waves in the FN kinetics studied in Ref. [26|. The kinetics is defined by the 
equations dtg = DVg — k r g(g — a)(g — 1) — k r g and dtr = (g — r)/r. We used k g = 5.2 and the 
other parameters as defined in Ref. pq|: D = 1, k r = 1.5, a = 0.05, and r = 5. 
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